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PREFACE TO THE SECOND EDITION 




This second edition of the IDSP documentation, in addition to a number of 
minor corrections and updates to the original 1982 docunentation , contains 
substantial new material. Since the Fall of 1982 five new operators have been 
added and are included in this edition. 

RECPOL provides for rectangular to polar, and polar to rectangular rotations, 
ROTATE provides for coordinate rotations in conjuction with, or independent 
of, the EIG operator. MEM implements the Maximum Entropy Method for computing 
power spectra. And FILOPT and WTS aid in producing an optimal digital filter 
design in conjuction with FILDES. 

Section 5.1 has been expanded to contain an example of the use of IDSP for the 
analysis of waves. This example demonstrates the use of operator EIG to 
perform an eigenfunction analysis on the components of a vector time series, 
Fourier transformed to the frequency domain, which reveals the polarization 
characteristics of the data as a function of frequency. In the process it was 
desirable to digitally filter the data and the example Includes information on 
how to optimize, using FILOPT, a two band digital filter designed with the 
Remez exchange algorithm. 

Section 5.2 has been expanded to discuss the use of FILOPT and WTS in the 
design of optimal digital filters and several examples are presented. See 
particularly Fig. 5.2-21 for an overview of the total filter design procedure. 



ABSTRACT 


The Interactive Digital Signal Processor (IDSP) is implemented on a VAX 

11/780 under VMS. It consists of a set of time series analysis "Operators" 
each of which operates on an input file to produce an output file ; the 
operators can be executed in any order that makes sense and recursively, if 
desired. The operators are the various algoritluns that have been used in 
digital time series analysis work over the years. In addition, there is 

provision for user written operators to be easily interfaced to the system. 
The system can be operated both interactively and in batch mode. 

In IDSP a file can consist of up to n (currently n=8) simultaneous time 
series. Thus storage for a file can be subdivided such that it is used, for 

example, entirely for one long single time series or for as many as n shorter 

time series, such as the components of a vector. An operator always operates 
simultaneously on all of the time series in a file. 

IDSP currently includes over thirty standard operators that range from Fourier 
transform operations ( FFT,FFTIN, WINDOW, SPECT) , design and application of 

digital filters (FILDES,FIL0PT, FILTER, WTS) , eigenvalue analysis (EIG), to 

operators that provide graphical output (GRAPH, GRAFCK) , allow batch operation 
(REDO), editing ( CONCAT, EDHIST, EDIT, INTERP, SUBSET, SUBSER) and display 
information (SHOW, CMDHIS) . The complete set of standard operators is listed 
below. 

AVER, CMDHIS, CONCAT, COPY, DCL, DSTAT, EDHIST, EDIT, EIG, FILDES, FILOPT, 
FILTER, FFT, FFTIN, FLOP, GRAFCK, GRAPH, INTERP, MEM, MNFLD, NORM, RECPOL, 
REDO, ROTATE, SETUP, SHOW, SPECT, SUBSER, SUBSET, TRACE, WINDOW, WTS, STOP. 

IDSP is being used extensively to process data sets obtained from scientific 
experiments onboard spacecraft such as Dynamics Explorer, ISEE, IMP and 
Voyager. In addition IDSP provides an excellent teaching tool for 

demonstrating the application of the various time series operators to 
artlfically-generated signals. 

IDSP is available from the Computer Software Management and Information Center 
(COSMIC), 112 Barrow Hali, University of Georgia, Athens, Georgia 30602, 
Program Number GSC-12862. 
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time. 
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(odd function) about the NyqulsX frequency . 
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Fig. 5.1-11 ' i 

Eigenfunction properties of fluctuations shown in Fig. 5.1-9 over a . - t}‘ . a 
range of frequency (4 X 10**-5 to 4 X 10**-4 Hz) that Includes the p*;. -ve 
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trace spectrum is repeated from Fig. 5.1-10 for reference. Also shown, in 
panels 2 through 4, respectively, are results from the application of EIG: 
degree of polarization, cosine of angle between % and It, and wave elliptlcity 
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Raw time series consisting of 9.6 sec averages for the portion of day 250 in 

which higher frequency fluctuations were seen also in Fig. 5.1-9 (1200-1824 
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magnitude (B„) is significantly lower than that of those in the vector 
components . 
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The transfer function (frequency response) of the high pass filter designed 
with operator FILDES to remove the low frequency oscillations from the raw 
data. The function was obtained and plotted through the successive operators 
WINDOW FILDES, FFT WINDOW, and GRAPH FFTWP. 

Fig. 5.1-15 Page 5.1-37 

Passband (PB) and Stopband (SB) ripple amplitrde as a function of relative 
bandweight WTX(1)/WTX(2) for a 125-coefficient high pass filter with the given 
bandedges. The vertical dashed line delineates ripple characteristics 
corresponding to relative weight value of 182 used for filter shown in Fig. 
5.1-14 (see text and Fig. 5.1-16). 
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0.0324 and 0.175 Hz. 
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The output of operator FILOPT showing 201og(1 t- delta.), and 201og(delta2) 
plotted against the weight ratio, delta ./delta-* From these plots we 
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Shows the results of using operator FILTER to filter the original time series 
shown In Fig. 5.2-2 with the 5 band filter shown In Fig. 5.2-10. 

Fig. 5.2-12 Page 5.2-19 
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Observation time of 0.5 sec. 
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Plot of a 10 Hz time series after being windowed by a Hamnlng window, note how 
the ends are tapered to zero. 
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Fig. 5.3-9 Page 5.3-9 
Transform of a 10 Hz time series contdining 1 unit of noise using a Hamming 
window and no zero padding. 


Fig. 5.3-10 Page 5.3-9 

Transform of a 10 Hz time series containing 1 unit of noise using a Hamming 
window zero padded out to 200 points, note that the leakage level is about the 
same as shown in Fig 5.3-9. 


I 


t 


i 

I I 

' 1 

i 

I 

I. 

I 


\ 1 


xii 


t 


list of figures 5/29/84 



1.0) INTRODUCTION TO IDSP 


Many, but of course not all, of the phenomena that we want to measure in the 
real world are of an analog nature. Often it is desirable, perhaps essential, 
that we be able to process these measurements on a digital device in order to 
be able to extract the desired information from the raw data. The process of 
mapping the s ■'■'og phenomena to the digital world is called sampling. When a 
phenomenon if ^mpled as a function of time, what results is a digital time 
series . 

The design of the instrumentation to perform this sampling is, of course, 
critical. How often the instrument samples the phenomenon (to avoid 
aliasing), the time between samples (is it constant?), the number of bits of 
digitization (precision), and the transfer function of the instrument all have 
a large affect on the ability of the researcher to analyze the data properly. 

Fig. 1.0-1 is a schematic of this process, i.e., an analog phenomenon is 
sampled as a function of time and digitized, and v^iat results is a digital 
time series. 

Often the phenomenon of interest requires that multiple parameters of the 
phenomenon be sampled similtaneously. As an example: if the phenomenon of 
interest is a vector quantity, then it is required that all of the components 
of the vector be sampled. This, in turn, will result in there being not one 
but ”n" digital time series generated. 

The processing of these digital time series to extract information usually is 
done with a number of fairly standard operations that are applied to the data 
serially, but, of course, not always in the same order. It may also be very 
desirable to "experiment" with various operations to explore their effect on 
the results. 

Some of the more common operations that come to mind are the editing of data 
to remove "bad" points, interpolation of the data during data drop-outs, 
digital filtering of data to remove unwanted signals, and discrete Fourier 
transformation of data so that it can be viewed in the frequency domain, to 
mention but a few. 

It was the above mentioned concepts that dictated the design of the 
Interactive Digital Signal Processor software. 

IDSP was designed and implemented to be extremely interactive and easy to use. 
Essentially, it consists of a set of "operators" each of which operates on an 
input file to produce an output file ; the operators can be executed in any 
order that makes sense with respect to the task to be accomplished and 
recursively, if desired. The operators themselves are simply the various 
algorithms that have been used in digital time series analysis work over the 
years. Some specific examples are: an operator FILTER to filter a time 
series, an operator FFT to perform a discrete Fourier Transform on a time 
series, an operator INTERP to interpolate across gaps in a time series etc. 
In addition, there is provision for user written operators to be easily 
interfaced to the system. 
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In order for an operator to process an input file to produce an output file, 
it is necessary that the input file have a "name". In IDSP the name of the 
input file is simply the name of the last operator to process this file. 
Thus, if we have just interpolated a file, the name of this file is now 
INTERP. If we now filter INTERP, its name becomes FILTER. In this fashion a 
file progresses through the various operations. Note that when the data are 
first introduced into the system the name of the file is SETUP. 

The only exception to this rule of the output file being named exactly for the 
operator that just processed it is the case where an operator produces more 
than a single output file. An example is the spectral density matrix 
operation SPECT which, for technical reasons, produces three output files 
SPECTD and SPECTOFF containing the diagonal elements of the matrix and off- 
diagonal elements, respectively and COHPH which contains the coherence and 
phase information arising from the off-diagonal terms. 

So far we have treated the terms "time series" and files as if they are the 
sane thing; however, there is a distinction. In IDSP a file can consist of up 
to n (currently n=8) simultaneous time series, as long as each tuple of the 
file has the same time tag (see Fig. 1.0-2). Thus storage for a file can be 
subdivided such that it is used, for example, entirely for one long single 
time series or for as many as n shorter time series (total space not to exceed 
500,000 real points, operations involving FFT will be smaller due to working 
arrays). An operator always operates simultaneously and identically on all of 
the time series in a file regardless of whether there is one, two, ..., or n. 
This distinction is an important aspect of IDSP when one considers that many 
time series consist of multiple components, e.g., vector components. 

To allow the user to keep track of what operations have been done to a file, 
an integral part of each file is hibcory information that explicitly records 
the history of operations performed. The history information is very useful 
because the processing of a file can involve many operators in a nontrivial 
application, and it is quite easy to forget just vrtiat operations have been 
applied. Any time the user wishes to find out this history it is only 
necessary to type in SHOW "file name" and the history of the file to this 
point will be displayed on the terminal. Also when a file is graphed with the 
Operator GRAPH the history information is displayed along with the graph. 

IDSP processes the time series according to two basic concepts: SPAN and 
INTERVAL (see Fig. 1.0-3). The Interval is defined as the basic time segment 
to be analyzed~the segment of data that is to be filtered or Fourier 
transformed, etc. by an operator. The Span is defined to be the total time 
under analysis — composed of one or more contiguous Intervals . 

IDSP is designed to operate on a DEC VAX 11/780 and thus the notion of 
DIRECTORIES needs to be understood. Since a disk on the VAX can contain files 
belonging to many different users, each disk has a set of files called 
directories, i.e., a catalog of the files on that disk that belong to a 
particular user. When the user "Signs On" the VAX the user is automatically 
put into his/her default directory. IDSP uses one sub-directory of ^this 
default directory: [Userid .IDSP] . The sub-directory [Userid. IDSP] is use^' to 
store input parameters required by some of the operators, e.g., in order to 
design a filter there is an operator called FILDES which requires that certain 
design parameters be supplied; these are stored in sub-directory 
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[Userid. IDSP], The same sort of thing is true of operator SETUP. 

[Userid. IDSP] is also used to store the executable version of IDSP and all of 
the vfor lei ng files of I^)SP ut.d it is out of this sub-directory that execution 
takes pic e. The execution of an operator, as previously explained, generates 
an output data set. Only the latest version of each data set is retained by 
IDSP, e.g., multiple execution of the same operator over-writes previous 
versions of the data set unless the COPY operator is used to save the previous 
version. 


IDSP includes the following Operators; 

1. AVER— average every n points 

2. CMDHIS — give listing of commands issued during current session 

3. CONCAT — concatenate two different datasets 

4. COPY — copy one dataset to another 

5. DCL — allow user to enter VAX DCL commands and return to IDSP 
Decimation — reduce data as option when filtering with operator FILTER 

6. DSTAT — compute dataset statistics 

7. EDHIST— allow user to edit history of dataset 

8. EDIT— display series and interactively edit on a HP2648A Graphics terminal. 

9. EIG — rotate spectral matrices to eigenvalue coordinates 

10. FILDES — design filter via Remez Exchange method 

11. FILOPT — aids in producing an optimal filter design 

12. FILTER— filter data 

13. FFT — discrete Fast Fourier T1"ansform 

14. FFTIN — Inverse discrete Fourier Transform 

15. FLOP— interactively select printout device during execution 

16. GRAFCK— check to see if sub process executing graph Job has finished 

17. GRAPH — plot each series in the dataset 

18. INTERP — interpolation 

19. MEN — Maximum Entropy method of computing power spectra 

20. MNFLD — transform a vector to mean field coordinates 

21. NORM — normalize data 

22. RECPOL — rectangular to polar; polar to rectangular rotations 

23. REDO — repeat command sequence in batch Jobs 

24. ROTATE — arbitrary coordinate rotation 

25. SETUP — call application interface INPUT, place data into proper 
configuration for analysis 

26. SHOW — display history and optionally data of given dataset 

27. SPECT — form spectral matrices 

28. SUBSER — extract specified series from dataset 

29. SUBSET — create a subset of a given dataset 

30. TRACE — form a trace of specified dataset 

31. WINDOW — choose data window, with option to pad with zeroes 

32. WTS — compute the number of digital filter weights 

33. STOP— termination procedure 

In addition to the above standard operators provision has been made for the 
user to write their own operators and interface them to IDSP (see Appendix B 
for details) . 
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2.0) AN EXAMPLE OF THE USE OF IDSP 


The following example illustrates how to use IDSP. On the Dynamics Explorer 
Spacecraft there is a magnetic field experiment that provides a vector 
measurement of the ambient magnetic field every 0.5 sec. TJie data from one of 
the components of the vector is shown in Fig. 2.0-1 as a function of time. 
This plot was obtained by using operator GRAPH on data set SETUP. 

As is evident there appears to be substantial periodic signal riding on top of 
the actual signal of interest. It was desired to remove this contamination 
from the signal of interest. First it was necessary to do a spectrum analysis 
on the data to determine the freque..cy(s) of the unwanted signals present. 
Fig 2.0-2 shows the results of using operators WINDOW, FFT, SPECT, and GRAPH 
on this time series, note that at 0.0324 and 0.175 Hz there are peaks in the 
spectrum resulting from the unwanted signals. 

We now use operator FILDES to design a finite impulse response (FIR) filter to 
remove both unwanted frequencies. In this case we use a single filter that 
has two stopbands and three passbands for a total of 5 bands. Operator FILDES 
will permit you to design a FIR filter with up to 10 bands total. 

The transfer function of this 5 band filter is shown in Fig. 2.0-3. This 
transfer function is obtained by using operator FFT on data set WINDOW which 
contains the windowed filter coefficients. The plot is obtained with operator 
GRAPH. Table 2.0-1 results from the execution of operator FILDES and presents 
the coefficients and other parameters of this 5 band filter. Fig. 2.0-4 shows 
the results of using operator FILTER to filter the original time series shown 
in Fig. 2.0-1 with this 5 band filter. 


The IDSP command sequence without required and optional operands is shown 
below: 


SETUP 

GRAPH SETUP 
WINDOW SETUP 
FFT WINDOW 
SPECT FFT 
GRAPH SPECTD 

FILDES 

WINDOW FILDES 
FFT WINDOW 


GRAPH FFTMP 
FILTER SETUP 
GRAPH FILTER 


input Dynamics Explorer data into IDSP 

plot the raw data 

window the raw data 

Fourier transform the windowed data 

generate the spectral matrices 

plot the transformed data to obtain the 

power spectrum 

design the 5 band digital filter 

window the filter coefficients 

Fourier transform the windowed 

coefficients to obtain the transfer 

function of the filter 

plot the filter transfer function 

digital filter the raw data 

plot the filtered data 
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Fig. 2.0-1 

On the Dynamics Explorer Spacecraft there is a magnetic field experiment that 
provides a vector measurement of the ambient magnetic field every 0.5 sec. 
The data from one of the components of the vector is plotted as a function of 
time. 
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INDEPENDENT VORIRBlE 

Fig. 2.0-2 

Shows the spectrum resulting from application of operators WINDOW, FFT, SPECT, 
and GRAPH on the time series from Fig. 2.0-1, note that at 0.032M and 0.175 Hz 
there are peaks in the spectrum resulting from the unwanted signals. 
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Fig. 2.0-3 


The transfer function of the 5 band filter is shown. This transfer function 
is obtained by using operator FFT on data set WINDOW which contains the filter 
coefficients. 
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Shows the results of using operator FILTER to filter the original time series 
shown in Fig. 2. 0-1 with the 5 band filter shown in Fig. 2.0-3. 
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3.0) OVERVIEW OF THE IDSP SYSTEM DESIGN 

1. The system must be interactive, yet capable of running in background mode. 
Basic functional operations are called operators. Each operator modifying a 
dataset writes results to a dataset identified with the operation performed. 
Other operators can then input those results by accessing the dataset 
associated with a previous operation. For routine operations, a command 
procedure can be written to batch process the data, using the operators in any 
meaningful order. 

2. The user chooses the operations appropriate for the analysis of his 
particular data, and is able to specify which dataset will be the inout for a 
given operator. If at any point in interactive mode the user feels a result is 
unsatisfactory, he can redo any operation using a dataset from an earlier 
operation as input. Or he may operate recursively by specifying the output 
from a given operator as input for that same operator. This structure allows 
the user to "interactively experiment" with the analysis of his data. 

3. Each dataset Includes records detailing the operations that have been 
performed on that dataset. This is included to help the user keep track of 
data set processing so that he does not misorder or duplicate operations (for 
example, inadvertantly filtering a dataset twice). 

4. The system is able to procr . up to n (currently n=8) different time series 
concurrently with up to NPTS p i 'ts each, as long as the total space required 
by the operation is less than a real array of length 500,000. NOTE: Because of 
the size of working arrays required, operations involving the FFT will be 
considerably smaller. For series involving more than 10,000 points, the user 
should check documentation under Design Considerations, APPENDIX C and also be 
aware that time of such long runs may be excessive. 

5. Any given operator imposes the same operation on all of the series 
currently being processed. Each operator assumes that every series being 
processed involves the lame time interval (see Fig. 1,0-2). 

6. The system processes the time domain according to two basic parameters: 
span, interval. The interval is defined as the basic time segment to be 
analyzed — the segmei.t of data which is to be filtered or Fourier transformed 
by a specified operation. The span is defined to be the total time under 
analysis — composed of one or more contiguous intervals (see Fig. 1.0-3). 

7. At the outset (in operator SETUP), each series is required to have valid 
data at the endpoints of the desired interval. If necessary, the system 
successively reduces the interval until this requirement is met. 

8. Cc.aparisons between series of different types where different operators are 
used can be accomplished by making separate runs, storing the results, and 
using this as input for a final run, A cross correlation operator (not 
currently available as a system operator) could h** >*::ed to analytically 
compare the two results. For meaningful comparison, each series must have the 
same number of points. A concatenation operator will be used when one desires 
comparisons between two different runs. This concatenation will be permitted 
only if the number of points in the two sets of series is Identical. This 
composite dataset can then be accessed as input for comparison operators. 
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9. The user can write his own operators ::ricl interface them into the IDSP 
system by adhering to the conventions described in Writing User Defined 
Operators, APPENDIX B. 

10. Each application accesses data by an input routine specxi^c tc the 
experimenter data of the application. This routine is called by the system, 
and is the interface between the system and the experimenter data. Because a 
different input routine is required if the application changes tne user will 
need to relink IDSP with this different input routine in this event. 


% 
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^. 0 ) HOW TO START USING THE IDSP 


The beginning user should do the following: 

1. Create a subdirectory called [usrid.IDSP] where the user will place the 
input parameter files. Use the following sequence: 

SET DEF [usrid] 

CRE/DIR [usrid.IDSP] 

SET DEF [usrid.IDSP] 

2. Create datasets in [usrid.IDSP] which are required by the SETUP command; 
FOR051.DAT (see operator SETUP, Section 6); FOR052.DAT (see Writing Input 
Routines* APPENDIX A), and any other data sets required by the i nput routine. 

3. ASSIGN experimenter data to F0R011 or other FORTRAN logical unit number as 
required by the specific input routine being linked (this step may be optional 
depending on the specific requirements of the INPUT routine being used). See 
Writing Input Routines, APPENDIX A. 

4. Activate the IDSP system by: @SYS$IDSP: IDSP linknames, where "linknames” 
are the names of the dataset(s) containing the object versions of special 
routines the user wishes to link into the system, e.g. user written operators. 
The name of the dataset containing the desired applications interface INPUT 
routine is always required; datasets containing user defined operators are 
optional (see Writing User Defined Operators, APPENDIX B). Multiple datasets 
must be separated by commas. The user should indicate the device name where 
appropriate to avoid a search of the wrong device and subsequent abort of the 
link. If the user has already linked his desired routines in a previous run, 
he may skip the link step by specifying "NOLINK" for "linknames". WARNING: In 
order for the system to be properly initialized this execution step must be 
entered after every time the user logs on. (Otherwise lack cf proper system 
assignments will cause the program to abort.) 

Assuming the default device is DRA1:, IDSP will then do the appropriate link, 
make logical assignments, and set the default directory to DRA1 : [usrid.IDSP]. 
When it is ready for the user to begin using the system, it will display: 
ENTER COMMAND. Files created by IDSP will be stored into this default 
directory. 

5. Help in understanding the commands available may be obtained by typing 
»HELP», or by typing ’HELP command’, e. g.. ’HELP SETUP’. The first non-HELP 
command of any session must be either SETUP (place experimenter data into 
proper configuration for analysis, see Section 6) or FILDES (design a filter 
and store coeficients, see Section 6) unless datasets have been saved from a 
previous session. The session is terminated by the command STOP which executes 
termination procedures. 

ALL COMMANDS MUST BE TYPED IN UPPERCASE. 


Commands have the following form: CMD REQ OPT where CMD is the command name, 
REQ are the required parameters (separated by spaces), and OPT are the 
optional parameters (separated by spaces). For any given command, the user 
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will be prompted if required parameters are omitted. Unless optional 
parameters are entered on the first line, their default values will be used. 


6. An example of a reasonable sequence of commands follows: 

SETUP (obtain interval for analysis) 

INTERP SETUP (interpolate missing points) 

GRAPH INTERP (Versatec plot of data) 

FFT INTERP (do Fourier transform) 

SPECT FFT (form spectral matrices) 

HELP SPECT (find out diagonal is stored in SPECTD) 
SHOW SPECTD (see operation history of diagonal) 
HELP GRAPH (recall graphing options) 

GRAPH SPECTD (see graph of spectral terms) 

STOP (end the IDSP session). 
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5.0) APPLICATION NOTES 


5.1) FOURIEH TRANSFORMS. SPECTRAL DENSITY MATRIX, THE EIGENVECTOR SYSTEM AND 
MAXMIUM ENTROPY METHOD (OPERATORS FFT, FFTIN, SPECT, EIG, MEM) 

5.1.1) OPERATOR FFT 

We start by considering three related concepts of Fourier analysis, t: 
Fourier Series, the Continuous Fourier Transform (Integral) and the Discrete 
Fourier Transform. 

It can be shovm that under rather liberal conditions (see Lanczos, 1956) an 
entirely "unp'^edictable" function f(t), generally normalized to the range ±n, 
can be represented, (to any arbitrarily high degree of accuracy), by a sum of 
components which are periodic (sines and cosines). This sum or series is 
called a Fourier Series . The representation of f(t) as a Fourier series 
demands strict periodicity in the time domain. That is, when f(t) is 
represented by a Fourier Series it is assumed to repeac with period 2ir outside 
of the fundamental range ±». The function f(t) may be a truly periodic 
function, or may exist in the finite interval only, and we force 
periodicity on it in order to make the Fourier Series applicable for its 
representation. In the latter case periodicity is employed as a mathematical 
artifice . 

Fourier found that decomposition of arbitrary functions into harmonic 
components remains possible even if the realm of the function f(t) extends 
beyond ±ir to ±® (see Lanczos, 1966). This i*» called the Continuous Fourier 
Transform (see Eq, 5.1-1) which is neither periodic in the time or frequency 
domain due to the infinite limits. This new function, X(f), does not resemble 
the original function in any direct way but is merely associated with it 
somewhat as the logarithm of a number is associated with the original number. 
For the purpose of time series analysis this process maps the function from 
the time domain to the frequency domain. 

If the data we have to work with are digitized observations taken at 
equidistant time intervals, AT, we employ the methods of Fourier IVansforms 
but adapted to finite summation ( Discrete Fourier Transform ) instead of 
Integration (Continuous Fourier Transform). The Discrete Fourier Transform is 
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a special case of the Continuous Fourier Transform where it is assuned that 
the N samples of the original function f(t) are one period of a periodic 
waveform. 

The Discrete Fourier transform, which is implemented as operator FFT, is of 
interest because it, under certain conditions to be discussed, approximates 
the continuous transform and thus allows us to Fourier transform discrete 
(digitized) data. The Discrete Fourier Transform implies periodicity in the 
time domain which results in periodicity in the frequency domain. 

The relationship between the discrete and continuous Fourier transform is as 
follows : 

The Continuous Fourier transform is defined as: 

+OD 

X(f) = 1/(2ir)/x(t)e"^^*^’^dt +® ; i=(-1)°*^ (5.1-1) 


The equivalent Discrete Fourier transform is defined as: 


N-1 

X(J)= j,k=0, 1,2 N-1 (5.1-2) 

k=0 


where N = the total nunber of data points in the original time series. 

The validity of this relationship is a function of the particular waveform 
being analyzed. 

According to Brigham, 197^» there are five cases of time domain functions to 
consider : 

Case 1 Band-limited periodic waveform, truncation interval (rectangular 
window) equal to one (or multiple) perlod(s). 
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This example represents the only class of waveforms for which the discrete and 
continuous Fourier transform are exactly the same within a scaling constant. 


Equivalence of the two tranforms requires; 

(i) the time series of interest, x(t), must be periodic over T (see Fig. 
5.1-1); 

(ii) x(t) must be band-limited; 

(iii) The sampling rate, f^, must be at least two times the largest 
frequency component of x(t); It is defined as f^=1/(2AT), where AT is a 
constant and is the time between samples; and 

(iv) The truncation (rectangular window) must be non-zero over exactly one 

(or multiple) period(s) of x(t). This also implies that the time series 
x(t) should be’ a stationary series. This means its statistical properties 
should not change with time during the period of time spanned by the 

series (see Jenkins, I968). 

Case 2 Band-limited periodic waveforms, truncation interval not equal to 
integer period(s). 

If a periodic, band-limited function is sampled and truncated to consist of 
other than an integer multiple of the period, the discrete and continuous 
fourier transforms will differ. The effect of truncation at other than a 
multiple of the period is to create a periodic function with sharp 
discontinuities (see Fig. 5.3-7 under the discussion of zero padding). These 

sharp discontinuities in the time domain result in additional frequency 

components in the frequency domain. This effect is termed leakage . Windowing 
with other than a rectangular window (see Section 5.3) can be employed to 

reduce this leakage. 

Case 3 Another class are functions which are of finite duration in the time 
domain. If x(t) is time-limited (for an exerjple see Matthaeus and (Soldstlein, 
1982), its Fourier transform cannot be band-limited and sampling must result 
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in aliasing. It is necessary to choose the sample interval AT such that 
aliasing is reduced to an acceptable range. For this class of functions, if N 
is chosen equal to the number of samples of the finite-length function, then 
the only error is that introduced by aliasing. Errors introduced by aliasing 
are reduced by choosing the sample interval, AT, sufficiently small and, in 
the limit (see Lanczos, 1956), the discrete Fourier transform vd.ll agree to 
v#ithin a constant vdth the continuous Fourier transform. 

Case 4 General periodic vjaveforms, truncation interval one (multiple) 
period( s) . 

Periodic waveforms not band-limited but truncated to an interval of one 
(multiple) period(s) will result in the discrete and continuous Fourier 

transforms oeing the same with thi only source of error being aliasing. If 
the truncation is not equal to an integer multiple of the period, then results 
are as described in Case 2. 

Case 5 General waveforms, not time-limited or band-limited. 

This important and commo'" class of functions are neither time or band-limited. 
Sampling thus results in an aliased frequency function and time domain 
truncation introduces rippling in the frequency domain. As this class of 

functions is often encountered vre would llKe to treat them as though they were 
either band-limited (Case 1) or time-limited (Case 3). The aliasing error can 
often be reduced to an acceptable level by decreasing aT, and the time domain 
truncation error can often be reduced by windowing with other than a 
rectangular window (see Section 5.3). 

Graphical interpretation of operator FFT 

For the following discussion see Eq. 5.1-2 and Fig. 5.1-1 which follows 

Bergland, 1969. Let (AT) be the time between samples in the time domain, so 

that the fundanental frequency is, f^sl/T, and f^ = Nf^. X(J), the discrete 

Fourier transform, is in general a complex series. The time series x(kAT) is 

assumed to be periodic in the time domain of period T. The Fourier 

coefficients X(Jf-) are periodic over f by definition of Eq. 5.1-2. Each J 
0 s 
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should be interpreted as a harmonic nunber and each k a sample period number. 
Note that x(NaT)=x(0), X(Nfp)=X(0), actual frequency = and actual time s 

kAT. 

When the series x(k) is composed of real nunbers, as it often is, the real 
part of X(j) is symmetric (even function) about the Nyquist frequency and the 
Imaginary part is antisymmetric (odd function) about the Nyquist frequency. 
Fig. 5.1-1 shows the relationship between x(kAT) and X(jfQ). This also can be 
seen in the example shown in Figs. 5.1-2 and 5.1-3, where we show the real and 
imaginary parts of the results of discrete Fourier traiiSforming x(t) = 10 Cos 
(2TT*15*k»AT) + 5 Sin (2ir*20»k»AT ) sampled at 120 times per sec. Note that the 
real part shown in Fig. 5.1-2 is symmetric about f^ (i.e., BIN 120/2 = 60) and 
the imaginary part shown in Fig. 5.1-3 is antisymmetric about f^. Fig. 5.1-4 
is the results of plotting the magnitude of this complex transform, which 
shows the expected peaks at BINS 15 and 20, respectively. Note: for 

convenience the range is thought of as 1 - to (N/2)fQ. 

Operator FFT computes the XCJf^) complex coefficients; only (N/2)+1 
coefficients are retained in the output data set (FFT) because of symmetry, 
where N is the number of points in the original time series (N=NPTS). 
Operator FFT also has the option of computing the magnitude and phase of the 
coefficients which are stored in data set FFTMP, if this option is executed. 

5.1.2) OPERATOR SPECT 

Often one wants to look at power spectra of a vector quantity Bj^(i=1 ,2, 3). 
for this it is convenient to look at the matrix G = G(jf(\) * 

<Bi ( jfjj) ,Bj( jfQ>>, which is an estimate of the Fourier transform of the 
2-time correlation matrix <Bj^(t) ,Bj(t+r)>. For any vector quantity, the 

matrix G can be written in terms of the Power Spectral Densities (PSD), 
cospectra and quadrature spectra associated with the three ve^'tor component 
directions (x,y,z) as follows (see Otnes and Enochson, 1972): 


Note: B^ = Complex ConJ. 
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(5.1-3) 


Where S is the power density and C is the cospectrun and Q is the quadrature 
spectrun. The sample cospectrum measures the covariance between in-phase 
components (i.e., between cosine components and sine components separately), 
and the sample quadrature spectrum measures the covariance between the 
out-of-phase components (i.e., between sine and cosine components, see 
Jenkins and Watts, 1968). 


Now C 


ij 




and 


‘ij 


-Q 


Ji 


by definition, so we may write 


G = 


( S, ) 
(^xy^iQxy^ 

^^xz^'«xz> 


(Cxy-'Qxy> 


(C +iQ ) 
yz yz 


-V“Vx> 


(5.1-4) 


or 
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Thus the PSD matrix is composed of a real, symmetric part and an imaginary, 
skew-symmetric part, and hence the matrix is hermitian. 


We can also use the off-diagonal terms of the power density matrix to compute 
coherence and phase lag. The coherence function is defined as 


^ij = (5.1.6) 

and indicates whether the amplitude of the component at a particular frequency 
in one aeries is* associated with a large or small amplitude at the same 
frequency in the other series. 
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and phase angle or phase lag in degrees is defined as 

♦jj = (180/ir)tan""' (Q.y'C.j) (5.1-7) 

and indicates whether the frequency components in one series lead or lag the 
components at the same frequency in the other series. 

Operator SPEC! computes the PSD matrix by using ^ to the first four (*1) 
complex series in the file generated by Operator FFT and in turn generates all 
possible permutations of the power and cross spectra using equation 5.1-8 
below; note that i is the index for the individual elements of each time 
series, and j is the series index. 

P(i,J,k) = XNORM»X(i,j)*«X(i,k) (5.1-8) 

Note: X(i,j) = Complex Conj., 




where XNORM = 2/( (2*NPTS-2)**2) , thus folding the negative power and adding it 
to the positive. For i=1, XNORM = 1/( (2«NPTS-2)»*2) . 


Example 


For the case of three (3) time series, each containing 11 points (NPTS=11), 
these matrices are generated in the following fashion. Note that the time 
series, after being Fourier transformed (Operator FFT), consists of 6 complex 
values of X(J). Note that for a given value of the time series index 1, Eq. 
5. 1-8 produces from the Fourier transformed i ’ elements of the three parallel 
series all possible permutations of auto and cross power spectra. 


P(1, 1, 1) 
P(1, 1,2) 
P(1.1.3) 
P(1,2,2) 
P(1,2,3) 
P(1,3,3) 
P(2,1,1> 


XN0RM»X(1, 1) 
XN0RM»X(1,1) 
XN0RM»X(1, 1) 
XN0RM»X(1,2) 
XN0RM»X(1,2) 
XN0RM»X(1,3) 
XN0RM»X(2, 1) 


•X(1,1) 

•X(1,2) 

•X(1,3) 

•X(1,2) 

•X(1,3) 

•X(1,3) 

•X(2,1) 
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P(5,3,3) 

P(6,1,1) 

P(6,1,2) 

P(6,1,3) 

P(6,2,2) 

P(6,2,3) 

P(6,3,3) 


XNORM«X(5,3) 

XN0RM»X(6,1) 

XN0RM»X(6,1) 

XN0RM»X(6,1) 

XN0RM»X(6,2) 

XN0RM»X(6,2) 

XNORM*X(6,3) 


•X(5.3) 

•X(6,1) 

•X(6,2) 

•X(6,3) 

»X(6,2) 

•X(6,3) 

•X(6,3) 


For the above calculations Operator SPECT would generate the set of PSD 
matrices shown below. Note that these Matrices are hermitian so we only 
compute the upper right portion: 


i=6 


P(6,1,1) P(6,1,2) P(6,1,3) 
P(6,2,2) P(6,2,3) 


P(6,3,3) 


is5 


P(5,1,1) P(5,1,2) P(5,1,3) 

P(5,2,2) P(5,2,3) 
P(5,3,3) 


i=4 


P(«,1,1) P(4,1,2) P(‘»,1,3) 
P(4,2,2) P(4,2,3) 
PC'*, 3,3) 


i=3 


P(3,1,1) P(3,1,2) P(3,1,3) 
P(3,2,2) P(3,2,3) 
P(3, 3,3) 


i:2 


P(2,1,1) P(2,1,2) P(2,1,3) 
P(2,2,2) P(2,2,3) 
P(2, 3,3) 


is1 


P(1,1,1) P(1,1,2) P(1,1,3) 
P(1,2,2) P(1,2,3) 
P(1,3,3) 
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Operator SPE^T partitions the PSD matrices into data set SPECTD containing the 
diagonal components which are real and provide information about the power at 
each frequency in each separate time series and data set SPECTOFF containing 
the off-diagonal elements which are complex and contain information about the 
phase and coherence between time series. In the above example P(i,1,1)» 
P(i,2,2), P(i,3»3> vd.ll be contained in SPECTD while the off-diagonal terms, 
P(i,1,2), P(i,1,3), P(i,2,3). v»ill be contained in SPECTOFF. Operator SPECT 
uses equations 5.1-6 and 5.1-7 to compute the coherence and phase lag. These 
values are stored in complex dat.3 set COHPH. The real part of COHPH contains 
the coherence and the imaginary part contains the phase lag. See Fig. 5.1-5 
for a schematic of Operator SPECT's input and output data sets. 

5.1.3) STATISTICAL STABILITY 

Part of Operator SPECf is an option to allow PSD matrices to be averaged 
together to improve statistical stability of estimates being made. 

X(j) is a measure of the power in a frequency range bounded, for the 
estimate, by the frequencies (j±1/2)fjyM and centered at jf^M where M is the 
number of spectral estimates to be calculated. 

Associated with the computation of PSD is the degree of freedom u, which for 
moderate to large values is given by; 

u = 2N/M (5.1-9) 

where N, as before, is the number of data points (NPTS) in the given time 
series. The quantity u is related to the confidence interval for the 

estimates (see Sentman, 1974). for u > 4, confidence limits can be computed 

with good accuracy using the confidence factor k , where 

c 

k^ s exp(2.3b/10(u-1)°*^) (5.1-10) 

and for 98S confidence, b s 29; for 96t, b s 25; for 901, b s 20; and for 80S, 
b 3 16, Using k^ one can calculate the error bar limits from 

upper limit s (PSD estimate) (k^)®*^ 

c 
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lower limit r (PSD estimate) /(k^)^‘^ 

For example, if N = 1500 samples and M = 60 spectral estimates, giving 

u = 2N/M s 50 

then k = exp(0. 0329b) and, for 90X confidence, (k.)'^*^ = 1.390 and 1/(k^)^*^ 
c c c 

= 0 . 719 . The PSD estimate would thus be multiplied by the factors 1.390 and 
0.719 to obtain error bar upper and lower limits, respectively. 

Note the relationship between data set SPECTD, when not averaged, and the real 
part of data set FFTMP (see Section 6,FFT Operator); SPECTD = 
2»(FFTMP/(2*J-1))»»2, where J = the number of points in the FFT including zero 
frequency. 

5.1.'!) OPERATOR EIG, EIGENVECTOR SYSTEM 

Users of spectral analysis techniques are often interested in applying them to 
the study of wave phenomena. In such cases it is usually necessary to 
investigate the fluctuation characteristics of the individual components of a 
vector, not only relative to a physically-defined set of coordinates but also 
relative to one which is defined by the directional properties of the 
fluctuations themselves. For the case of an ideal elliptically-polarized 
plane wave, for example, there is no fluctuation perpendicular to the wave 
front, and a maximum level of fluctuation along one direction parallel to the 
wave front. There is thus a preferred coordinate system which most perfectly 
reveals the nature of such a wave. 

Such an analysis is Implemented using an eigenvalue and eigenvector 
calculation. In such a calculation, the characteristic variance ellipse (see 
Fig. 5.1-6) is determined, to the extent possible; the principal axes of the 
ellipse are the directions of maximum, intermediate and minimum fluctuation, 
respectively. To obtain eigenvalues and eigenvectors one needs only to 
diagonalize the real part of the PSD matrix, i.e., the real part of Eqn. 
\5.1-5) For this purpose, IDSP employs the subroutine EIGRS, which is an IMSL 
(International Mathematical and Statistical Libraries Inc., see Edition 8, 
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Vol. 2, Chapter E) FORTRAN routine which computes and returns the eigenvalues 
and associated eigenvectors of a real symmetric matrix. The eigenvalues are 
the variances in the three principal axis directions (see Fig. 5.1-6). EIGRS 
does not order the eigenvalues and eigenvectors according to maximum, 
intermediate and mimimum, however; this is accomplished in subroutine DIAG. 
The eigenvectors Vj are labeled with the indices JMAX, JMID, and JMIN, 
respectively, in the ordering process. When used in conjuction with spectral 
analysis, this process is carried out separately and independently for each 
spectral estimate (each frequency bin). Thus the directions of interest can 
be different at each frequency and Fig. 5.1-6 is for a particular frequency, f. 

For the computation of parameters which characterize the polarization 
properties of the fluctuations (see Fig. 5.1-7), we transform to a coordinate 
system defined by the eigenvector system. This can be done in one of two 
ways, either using the standard eigenvector definition in terms of the 
principal axis system or using a special definition appropriate for data in 
"mean field coordinates'- , which are defined in tne description of operator 
M^FLD (see Section 6). These systems are further described below. 

1) Standard Eigenvector System - Under this option (the default mode), tne 
standard eigenvector definition developed above is used. 

From this definition a matrix R can be derived for transformation of data from 
the input coordinate system to the eigenvector system: 


R 


V(1,JMAX) 
V(1, JMID) 
V(1, JMIN) 


V(2, JMAX) 
V(2, JMID) 
V(2, JMIN) 


V(3, JMAX) 
VO.JMID) 
V(3,JMIN) 


(5.1-11) 


Where V(i,JMIN), i = 1,2,3, is the unit vector in the minimum variance 

direction. 

This matrix is then used to transform both the real and imaginary parts of the 
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PSD matrix data sets (SPECTD and SPECTOFF). Polarization parameters are 
computed according to the formulation described under option (2) oelow. 

2) Mean Field Eigenvector System - Use of this option is valid only if the data 
analyzed have been transformed to the Mean Field Coordinate System by use of 
operator MNFLD (see Section 6) prior to application of EIG. This will have 
set the switch MFIELD = 1, which automatically activates the appropriate 
eigenvector derivation procedure. 

In this case, the ordering of eigenvalues and eigenvectors in DIAG will lead 
to the definition 


= V(1,JMIN) 

= V(2,JMIN) (5.1-12) 

kj = V(3, JMIM) . 

The unit vector k gives the wave normal vector direction in the case of a wave 
analysis, and is here taken as the direction of minimum fluctuation, l.e., the 
direction given by the eigenvector assiclated with the minimum eigenvalue, and 
is one of the set of eigenvector coordinates. 

For data in mean field coordinates, the matrix R can be derived which 
transforms the real and imaginary parts of the PSD matrix data sets (SPECTD 
and SPECTOFF) to (mean field) eigenvector coordinates. In this case, R is 
expressed in terms of the components (k^,k 2 *k^) of K only; 




-^ 2*^3 


VkfTkl" V^kJ ♦ 




R = 


2 _ 

1 * '^2 
k. 


-k < 




V^k J ♦ k| 


‘1 


(5. ;-i3) 
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This nratrix is then used to transform both the real and Imaginary parts of the 
PSD matrix data sets (SPECTD and SPECTOFF) . To compute polarization 
parameters (assiming a plane wave) , one uses only the components of the 
resultant spectral density matrix corresponding to the plane perpendicular to 
the wave vector direction. Thus the analysis reduces to that of a 2 x 2 
matrix. 


For the polarization analysis we have followed the formulations of Fowler 
et al., 1967 and Rankin and Kurtz 1970. The two-dimensjcnal power spectral 
matrix is given oy 


J = 


— 



— 


J 

J 


R 

(R - il ) 

XX 

xy 


XX 

xy xy 

J 

J 


(R 

- il ) R 


yy 


yx 

yx yy _ 


(S.l-I'*) 


where R and I signify real and Imaginary parts, respectively. Now |Jl 

R R - C(R - il ) (R - il )] or, since 

XX yy yx yx xy xy 


then 


R r R and I s - I , 
yx xy yx xy 


IJI = R R - R^ - . 

' ' XX yy xy xy 


(5.1-15) 


Computing the characteristic equation 


1 ? 1 /? 
Dr {(R + R ) - [(R + R )‘^ - t»| J| ] 

V V V V \/l/ ' ' 


XX yy 


XX yy 


from the characteristic root of the matrix J (see the derivation in Fowler et 
al., 1967), we can then define a matrix that represents the polarized part of 
the signal only: 


P r 


P„« P 


(J - D) J 

XX xy 

II 

XX xy 

P P 


J (J - D) 

yx yy 


yx yy 


(5.1-16) 


We then have that P r R •> D and P = R - D, and we can calculate the 

XX XX yy yy 

degree of polarization , DEGPOL, which is the ratio of polarized power to total 
power ill fluctuations, according to 
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where Tr 
coherency ; 


Tr|P| P + P,„, 

' ' _ XX yy 

Tr 1 Jl J n 

' ' XX yy 


(5.1-17) 


= trace. Further, we can compute, as previously shown, the 


J J 


xy 

yx 1 

J 

J J 

XX 

yy/ 


+ I 


sy xy 

1 T 

XX yy 


(5.1-18) 


The angle of polarization 6, defined as the angle between the X axis of the 
coordinate frame and the major axis of the polarization ellipse, is given by 


tan 2e 


J 

XX yy 


J - J 
XX yy 


where Re = real part. 
Thus 


THETA = 0 = ^ tan’^ 


J - J 
XX yy 


(5.1-19) 


For data in the eigenvector coordinate frame, THETA = 0. 

The ellipticity ELIP of the polarization ellipse is the ratio of the minor 
axis to the major axis of the ellipse. The parameter is computed as follows: 


sin 2s 


[(J + J )^ - 'll J| 

XX yy ' ' 
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where Im = imaginary part. Thus 


/ 


o ■> ^ -1 

6 = sin 


and 


-21 


xy 


l=<“xx - 


( 5 . 1 - 20 ) 


( 5 . 1 - 21 ) 


ELIP = tan |b|. 

The sign of b gives the sense of the polarization relative to the vector field 
being analyzed. This is shown, along with the corresponding phase angle 
value, in the Table 5.1-1. 


TABLE 5.1-1 



Le ft-handed 

Right-handed 


Polarization 

Polarization 

Phase lag 

90® 

270” 

R • 6 positive 

B negative 

6 positive 

R • 6 negative 

B positive 

B negative 


From thJs table we see that we can check results on the sense of polarization 
relative to the vector field at a given frequency for consistency between the 
computed phase angle at that frequency and the sign of B (dependent on the 
sign of k • S) . Since B is not one of the output parameters available in data 
set EIGPARM (see last paragraph of this Section), the sign of B has been 
assigned to DEGPOL, which is otherwise a positive definite quantity. Only the 
absolute value of B can be obtained from the quantity ELTP. The value of 
= BDOTK is an additional output quantity which depends for its sign on the 
sign of k. The sense of k (whether it is positive or negative) is completely 
arbitrary as determined by subroutine EIGRS. Thus it is necessary to adopt a 
convention for the sign of k. This is done by .leans of the optional parameter 
KCON. It allows the user of operator EIG to either force 1< always to be 
outward rather than inward-directed relative to the sun or to force Ic always 
to have a component in the +5 rather than the -6 direction (unless K ^ B) . A 
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more detailed description of KCON is given in Section 6.9. 

A second unit wave vector Ic is independently computed using the method 
described by Means, 1972. This technique derives the components of (< using 
only the imaginary part of the spectral density matrix. Good agreement 
between the results from the two methods has been taken as an indicator that 
the wave normal vector was well determined. The analysis output also includes 
the signal-to-noise ratio defined by Means, 1972 as 

SNR=Jll^. (5.1-22) 

•’zz 

where J ♦ J is the trace of the 2-ditnensional matrix defined by equation 
XX yy J 

(5.1-14) and J is the additional diagonal element of a 3 x 3 matrix that 
includes the l< direction. 

In studies of magnetic turbulence, a useful quantity (because it is an 

invariant of ideal MHD turbulence) is the magnetic helicity. The magnetic 

helicity is a measure of the topographical linkage of the magnetic field and 

is closely related to the polarization (see Moffatt, 1978). Techniques of 

determining the magnetic helicity can be found in Matthaeus and Goldstein, 

1982 and Matthaeus et al., 1982. If we denote by H (f) the reduced magnetic 

m 

helicity spectrun, then a parameterzation of bounded by *1 can be defined as 

SIGMA = fH (f)/tr|G| (5.1-23) 

m 

where f is frequency and G and are from Eq. 5.1-3 and H^(f) = 2Qy^(f)/f, 
where it is assuned that the spectrum has been reduced to the z axis (see 
Matthaeus and Goldstein, 1982 and Matthaeus et al., 1982 for more details). 
Note: Do not use the mean field option if you want to compute SIGMA- results 

will be meaningless . 

Data set EIGPARM contains scaler quantities DEGPOL, COH, ELIP, THETA, TRTOT, 
SNR, BDOTK, and SIGMA. EIG also creates data set EIGVEC which contains in 
series 1-3, the first three components of the eigenvector corresponding to the 
direction of minimum variance for that mode. Series 4-6 are the eigenvalues 
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(minimun, intermediate and maximun, respectively). Also created by operator 
EIG is data set EIGVEC1 which contains the three components of the eigenvector 
corresponding to the intermediate direction (series 1-3) and the three 
eigenvector components corresponding to the mc.ximum direction of variation 
(series 4-6). For a schematic of operator EIG and it data sets see Fig. 5.1-8. 

5.1.5) EXAMPLE: THE USE OF IDSP IN THE ANALYSIS OF WAVES IN DEEP SPACE DATA 

In this section we provide an example of the application of the IDSP to a 
specific scientific research task. In this case use of IDSP assisted in 
analyzing magnetic field measurements made by the Voyager 2 spacecraft 11 days 
after its closest approach to Saturn on August 27. 1981. At that time, low 
frequency waves were detected in the magnetic field, apparently in association 
with enhanced flows of charged particles from the direction of Saturn. Fig. 
5.1-9 shows the magnitude B of the total magnetic field as well as that of 
each of the three orthogonal vector components of the field in a heliographic 
coordinate system during a 24-hour period. The data plotted are 48-second 
averages of Initially 60 millisecond measurements. Large amplitude waves with 
a period of approximately 1.6 hours on average can be seen both in the 
magnitude and each of the components of the field. As noted on Fig. 5.1-9, 
higher frequency fluctuations of the field were detected also during a six 
hour period beginning at midday. They will be investigated later in this 
example. Ihe analysis of both types of waves was accomplished via the IDSP 
and will be illustrated. 

The data shown in Fig. 5.1-9 were processed using operators SETUP, INTERP and 
GRAPH. Using operator WINDOW, the effect of different types of windows on the 
data was examined. In the examples shown, the rectangular (default) window 

option was used. The data were spectrum-analyzed using the FFT and SPECT 

operators, with 3 points forming each spectral estimate in the latter 
operation (6 degrees of freedom). Fig. 5.1-10 shows part of the output of 
SPECT in the form of the power spectral density matrix diagonal terms plotted 
using GRAPH. They constitute the spectra of the field vector components and 
the field magnitude (data set SPECTD) . Shown along with the PSD of B for 
comparison is the trace of the PSD matrix (calculated using operator TRACE), 
which is a composite of the individual component spectra. The peaks of each 
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of the spectra, corresponding to the quasi-sinusoldal waves in the time 
domain, occur in a band centered at a frequency of 1.75 X 10 Hz. 

Ihe properties of the waves were studied further by means of the EIG operator. 
A subset of the results are shown in Fig. 5.1-11. Selected parameters were 
plotted from data set EIGPARM using operator GRAPH. Their values are shown 
for a restricted region of the frequency domain near the peak frequency, where 
in this case, a linear frequency scale has been chosen for display purposes. 
For reference, the corresponding segment of the PSD trace is plotted in the 
top panel. The most relevant quantities resulting from the eigenfunction 
analysis are shown below the trace: (1) Degree of polarization (DEGPOL), (2) 
cosine of the angle between the wave normal vector Ic and the magnetic field 
vector § (Coso = R • = BDOTK where o is known only modulo 180*-see 

discussion following Table 5.1-1/, and (3) the elllpticity of the polarization 
ellipse (ELIP). The definitions of these quantities were given in Eqns. 
5.1-17, 5.1-21 and the text following those equations. 


A vertical band one spectral estimate wide delineates the relevant region on 
each parameter curve in Fig. 5.1-11. Note that DEGPOL has been plotted with 
both positive and negative values. This is because this otherwise 
positive-definite quantity is multiplied by the sign of B (Eqn. 5.1-20). The 
wave is seen to be highly polarized ^DEGPOLp^gj^ s 0.95). An average value of 
BDOTK over the peak gives a = 112°; thus the wave is propagating at a large 
angle to the magnetic field direction. Positive values of DEGPOL at the wave 
frequency, together with negative values of BDOTK, suggests left hand (LH) 
polarization in the spacecraft frame of reference (see Table 5.1-1 and 
following text). However, as indicated following Table 5.1-1, in running EIG a 
convention forces Ic to have a particular sense of direction. The user must be 
careful that this direction is physically reasonable in terms of its 
implications for the specific phenomenon being studied. In this case, as 
would be usual for studies of waves in the solar wind, the analyst set KCON s 
1 (away from the sun). This also forced Ic to have a component toward Saturn, 
the source of the beam of ions exciting the waves. But, since the waves 

should have the same sense of direction as that of ion flow in the beam, (c 
should in fact be directed away from Saturn. Thus we are led to the 
Intepretation that in this case the wave was actually right hand (RH) 
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polarized (corresponding to a = 180* - 112* = 68®). Finally, an average value 
of FLIP : 0.7 is found at the peak, indicating that the wave is much closer to 
being circularly polarized (FLIP = 1.0) than linearly polarized (FLIP = 0). 

It is sometimes useful to replot the original data in a coordinate system 
defined by the set of oigenvectors corresponding to the directions of maximum, 
intermediate and minlmu.' oscillation in the wave. This is accomplished using 
the ROTATF operator on fie original (in this case INTFRP) data set, with the 
new coordinate system provided by the FIGVEC and FIGVFC1 data sets. The 
specification given in this case was: ROTATF INTFRP 1,2,3 0 NFRFQs16. The 
resulting data are shown in Fig. 5.1-12, with Xg being in the direction of 
maximum variation and Zg in minimum variation direction. 

As mentioned earlier, in Fig. 5.1-9 and 5.1-12, an interval of higher 
frequency fluctuations is seen superimposed on the 1.6 hour period waves. The 
msonetic field data for the 6.4 hour interval in which this second wave mode 
ocurred are plotted on an expanded time scale in Fig. 5.1-13. In this case 
higher time>resolution data (9.6 sec averages) are shown. Approximately 4 
cycles of the low frequency oscillation are still seen modulating the higher 
frequency fluctuations, which can be seen to be of much higher amplitude in 
the vector components than in the field magnitude. To facilitate the study of 
higher frequency variations, we removed the low frequency excursions by high 
pass filtering the data set shown in Fig. 5.1-13. Hie procedure for designing 
and applying digital filters is discussed and Illustrated in Section 5.2. The 
filter specifications that were entered into data set F0R058.DAT and used by 
operator FILOFS for this task are as follows: 


125,1,2,0 


0.0,0.002,0.02,0.5 


0 . 0 , 1.0 


182.0,1.0 


This is a filter with 125 coefficients, of the multiple passband/stopband 
type, with one stopband and one passband (for which the start ana stop 
normalized frequencies ar« (0.0, 0.002) and (0.02, 0.5), respectively). The 
grid density is set at 16 (input specification of 0 gives default value of 
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16) . Ideally there is complete attenuation in the stopband and no attenuation 
in the passband ( 0 . 0 , 1 . 0 ), and, finally, a relative weighting of 182 in the 
stopband to 1 in the passband is used. 

Fig. 5.1-1^ illustrates the frequency response of the filter generated for 
this task. Because the higher frequency oscillation to be studied was still 
at a relatively low value of normalized frequency (v> 0 .O 3 ), the passband 
response was designed to rise steeply to a value of 1 near zero normalized 
frequency. This requirement tightly constrained the values that could be 
given to most of the filter design parameters. The bandedge array values 
(EDGE(i), 1=1, H) and band weight functions (WTX( 1 ), WTX(2)) could be adjusted, 
however, within the given constraints, to permit to some degree the 
maximization of attenuation in the stopband and minimization of ripple 
amplitude in the passband. Because of the fact that the two factors work 
against one another for a given filter configuration as the weights are 
varied, reduction of passband ripple beyond a given point is limited by the 
onset of incomplete stopband attenuation. Conversely, stopband attenuation is 
limited by increasing distortion of the spectrum by the increase in passband 
ripple amplitude. Any such distortion is undesirable and may not be totally 
eradicable if the number of filter coefficients cannot be increased further 
and/or the bandwidths cannot be further adjusted, including the width of the 
transition band (see Fig. 5 . 2 - 1 ). 

A study was performed (using operator FILOPT, see Section 6) to provide 
insight into the degree of tradeoff possible in this example between stopband 
(band 1) attenuation and passband (band 2) ripple. While all other design 
parameters were held fixed, WTX(I) and WTX(2) were varied over values ranging 
from (1,100) to (1000,1), which, since the weights actually represent a ratio, 
is equivalent to setting WTX(2) s 1 and varying WTX(1) over the range 0.01 to 
1000. As described in Section 5 . 2 , for a given set of parameter values, 
FILDES generates a filter and lists its characteristics (as illustrated in 
Table. 5.2-1). Among the output characteristics tabulated are the amplitude 
deviations in each band, i.e., the delta^ and delta 2 values. Although the 
resulting values of each have been influenced by the full suite of input 
parameter values, and hence desired values may not be realized, the ratio is 
predetermined by the bandweights, i.e., by definition delta^/deltag = 

Page 5.1-20 


section 5.1 5/23/84 



WTX(1)/WTX(2) (see Rablner end (k)ld, 1975, PP 197.199). 


From the amplitude deviations, the passband and stopband ripple amplitudes In 
dB are computed (see Section 5.2 for definitions). 

In our tradeoff study we plotted the various passband and stopband ripple 
combinations output by FILDES is functions of the relative weight 
WTX(1 )/WTX(2) . We used the dB values since they are directly comparable to the 
desired performance specifications that enter as inputs into the filter 
planning process (see Section 5.2). The resulting curves are shown in Fig. 

5.1- 15. As can be seen, passband (PB) ripple is relatively flat over much of 
the region of steepest increase (negative) in stopband (SB) ripple. For 
larger relative weights 0160), however, PB ripple rises steeply, and the 
continued SB ripple growth becomes less steep. 

For our high pass filter, the stopband minimum frequency response (and hence 
the level of attenuation) is that computed at zero frequency (see Fig. 

5.1- 14). This value is also a function of relative bandweight, and, although 
influenced by the stopband ripple amplitude, its variation with changes in 
relative weight is different from that of the ripple amplitude. The minimum 
response (maximum attenuation) dependence on weight for this case is 
illustrated in Fig. 5.1-16. The curves shown in Fig. 5.1-15 and 16 can be 
obtained by using INTERP on the output from the FILOPT operator and graphing 
the results. 

Although there are two distinct minima in the curve, the second, occurring at 
a 'elative weight of 182, was found to provide the most complete attenuation 
of the low frequencies in the data. From Fig. 5.1-15 it is seen that for a 
weight of 182 the passband ripple has not yet grown large (0.13), but the 
stopband ripple is relatively large (-82), as desired. This was thus 
considered to be the optimum state for the given bandedge configuration. It 
is worth noting that prior to selection of bandedges (0.0,0.002,0.02,0.5), we 
investigated the configuration (0.0,0.001,0.01,0.5). For that case the ripple 
vs. relative weight diagram was similar in general form to that of Fig. 

5.1- 15, but differed considerably in detail. The stee Increase in PB ripple 
began at a low value of relative weight (4.5). The optimum values of SB and 
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PB ripple were -40 and 0.35, respectively, considerably less effective than 
those possible with the final bandedge configuration chosen. Thus it is also 
important to optimize bandedge values as well as band-weighting factors. It 
should be noted that a set of desired stopband and passband ripple amplitudes 
will end up being the achieved characteristics of a particular filter 
configuration only if the total set of design parameters can accomodate the 
realization of those values. 

Upon completion of the tradeoff and optimization procedure for our high pass 
filter design, operator FILTER was called to apply the resultant filter to the 
Voyager data. The filtered data set is plotted in Fig. 5.1-17. It is clear 
that the filter has successfully removed the modulation. Ihe remaining 
fluctuation is quite variable in amplitude, and the frequency is evidently 
variable as well. This may be the result of a superposition of several 
fluctuations of different a"'plitudes and frequencies. The corresponding 
spectral characteristics ' jhown in Fig. 5.1-18, irtiere 5 points were used to 
form each estimate (1C degrees of freedom). The rising trend of the PSD at 
low frequencies < 10 Hz in each case is the result of the filtering of the 

-3 

data. The weakly dominant frequencies in the fluctations (3.4 to 3.9 X 10 

Hz; period = 4.3 to 4.9 min) are denoted by the vertical dashed line; as can ■, 

be seen, this peak is only a single feature of a broad, highly structured • 

shoulder on the overall speocrum, as was anticipated from the appearance of 

the data 

* 

Selected results of the eigenfunction analysis are plotted in Fig. 5.1-19 

together with a , /rtion of the trace of the PSD matrix. As in the case of the 

1.6 hour period waves, the 5 min fluctuations are highly polarized (Maximum I 

DEGPOL s 0.8) and elliptical (Maximum ELIP = 0.75), polarization in the 

spacecraft reference frame is righthanded (RH), as indicated by +DP and 

-t-BDOTX, and the waves are propagating at a fairly large angle to § (average 

BDOTK :: cos a 0.57 do a s 55”). Note that in this case KCON = 0 was 

specified in running EIG, forcing Ic to have a component in the direction. 
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5.1.6) OPERATOR MEM (MAXIMUM ENTROPY METHOD ) 

The maximun entropy method (Operator MEM) for spectral estimates was 
introduced by Burg in 1967. Its purpose was to increase the spectral 
resolution when the length of the observed data record was short. 
Conventional spectrum analysis techniques usually window the autocorrelation, 
append zeros and then Fourier Transform to obtain the spectral estimate. 
Instead the MEM extrapolates the autocorrelation beyond the limited range of 
the data to derive the spectral estimates at the desired resolution with 
little or no interference from other frequencies. 

The MEM power spectrin, P(f), can be described by the following expression 


P(f) 


P At 
m 


(' 


1 = 


n.l 


mn 


) 


P^ = output power of the prediction error filter 

a 's = m+1 prediction error filter coefficients 

mn 

At = data sampling rate 

f s fVequency 


For details on how the filter coefficients are calculated see Anderson, 1974. 


The equation is solved by stepwise iterative increase of the matrix 
dimension from m-1 to m. The only remaining area of concern is choosing the 
number of filter coefficients, m. If ra is too small the power spectrum tends 
to be very smooth, thus eliminating any additional resolution; if m is too 
large, spurious information may be introduced into the spectrun. The value 
2N/in(2N) can be used as a default parameter. This was proposed by Berryman, 
1978 based on empirical evidence. Otherwise a trial and error approach is s 
tedious alternative. 
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fo= l/T ; fg=Nf ; T = NaT 


FIG. 5.1-1 

This drawing is after Bergland, 1969 and shows what happens when a time series 
x(t) is Fourier transformed. X(J) is, in general, a complex series. The time 
series xOc^T) is assumed to be periodic in the time domain of period T. The 
Fourier coefficients X(JfQ) are periodic over f . Each J should be 
interpreted as a harmonic number and each k a sample® period number. Actual 
frequency = jf^. Actual time s kAT. When x(k) series is composed of real 
numbers, as it often is, the real part of X(J) is symmetric (even function) 
about the Nyqulst frequency f-sf /2 and the imaginary part is antisymmetric 
(odd function) about the Nyquisc frequency. 
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Shows the real part of the results of Fourier transforming x(t) s 10 Cos 
(2«*15*lc*tT} * 5 Sin (2«*20*k*AT) sampled at 120 times per seo. Note that the 
real part is symmetric about the Nyqulst frequency, f^ (l.e., BIN 120/2 s 60). 
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j, 
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60.00 00.00 100.00 IjO.OO ^<40.00 

bins 


FIO. 5.1-3 

Shows the imagiRary part of the results of Fourier transforming x(t) • 10 Cos 
(2«*l5*k*AT) ♦ 5 Sin (2«*20*k*AT) sampled at 120 times per seo. Note that the 
Imaginary part is antisymmetrlo about the Nyqulst frequenoy, f^. 
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Fig. 5.1-4 

Shows the results of plotting the magnitude of the Fourier transform of x(t) = 
10 Cos (2n»15»k*AT) + 5 Sin (2ir«20»k»AT) sampled at 120 times per see, which 
shows the expected peaks at BINS 15 and 20, respectively and is symmetric 
about the Nyquist frequency. 



Fig. 5.1-5 

A schematic of the operators FFT and SPECT and their associated data sets. 
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Fig. 5.1-6 

In the study of wave phenomena it is common to investigate . using an 
eigenfunction analysis, the fluctuations characteristics of the individual 
components of the vector relative to the directional properties of the 
fluctuation themselves. Ih such calculation the eigenvalues determine the 
principal axes (oi*,ot*«ot*) of the characteristic variance (polarization) 
ellipse at each frequency estimate. The eigenvectors define the coordinate 
system corresponding to the directions of maximum (X), intermediate (Y) and 
minimum (Z) oscillation in the wave at each frequency estimate. 
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Fig. 5.1-7 
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Fig. 5.1-9 

Raw vector magnetic field time aeries from the Voyager 2 spacecraft consisting 
of 48-sec averages of field magnitude B„i and vector components Xjj, and Z„ 
in heliographic coordinates (ordinate units are nX s nanotestlas s 10**-5 
Gauss; abscissa units used are sample mmber rather than time). The data 
cover a 24-hour period (1981 /day 250) and Srow oscillation of field at two 
frequencies differing by a factor of 20. Note that at the lower frequency 
there are large amplitude variations in the field magnitude as well as in the 
components. 
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Fig. 5.1-10 

Power spectra for each time series shown in Fig. 5. 1-9i plus the trace (Tr) of 
the PSD matrix (upper left, with field magnitude B spectrum). Ordinate units 
are nTVHz and abscissa units are Hz. In each of the spectra a pronounced 
peak is seen centered on 1.75 X 10**-4 Hz (dashed vertical line) corresponding 
to the low frequency, large amplitude variation prominent in the time series. 
Tliis figure illustrates that there is relatively more power in directional 
fluctuations than in magnitude (field strength) fluctuations. 
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Eigenfunction properties of fluctuations shown in Fig. 5.1-9 over a restricted 
range of frequency (H X 10»»-5 to H X Hz) that includes the peak wave 

frequency (delineated by the vertical hatched band). In the top panel the 
trace spectrum is repeated from Fig. 5.1-10 for reference. Also shown. In 
panels 2 through 4, respectively, are results Qron the application of EIG: 
degree of polarization, cosine of angle between B and 1^, and wave elliptioity 
(see text for discussion). Page 5.1-33 
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Fig. 5.1-12 

Raw time series of Fig. 5.1-9 rotated into new coordinate system defined by 
the wave eigenvector set. Operator EIG generates sets of eigenvalues and 
eigenvectors for each spectral estimate. Those determining the coordinates of 
the data in this figure corresponded to the estimate containing the spectral 
peak denoted in Figs. 5.1-'!0 and 5.1-11. The Xg direction is that of the 
eigenvector associated with the largest eigenvalue, Z- is that associated with 
the smallest, and Yp is that associated with the intermediate value. Since 
the wave is nearly circularly polarized, the amplitudes of the oscillations in 
the Y- direction are nearly as large as in the Xg direction and essenti»M.< 
zero in the Zg direction. The magnitude of the field is not includ^.i" i.noc, it 
is invariant Aider rotation. 
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Raw time series consisting of 9*6 sec averages for the portion of day 250 in 
which higher frequency fluctuations were seen also in Fig. 5*1^ (1200-1824 
UT). Note that in this case the amplitude of variations in the magnetic field 
magnitude (Bj.) ^a significantly lower than that of those in the vector 
components* 
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Klg. 5.1-U 

The transfer function (frequency response) of the high pass filter designed 
with operator FILDES to remove the low frequency oscillations from the raw 
data. Ihe function was obtained and plotted through the successive operators 
WINDOW FILDES, FFT WINDOW, and GRAra FFIMP. 
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Fl«. 5.1-17 

Output from th# appllcutlon of th« high pass filter to the time “^ies show 
In Fig. 5.1-13. The low frequency modulation of the data has been removed 

successfully. 
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Fig. 5.1-19 

Eigenfunction properties of the fluctuations shown in Fig. 5.1-17 in the 
frequency band 0.001 to 0.01 Hz. Paraneters plotted are sane as those in Fig. 

5.1-11. Fk*equency of peak power in the fluctuations is delineated by the 
vertical hatched band. Conparison with Fig. 5.1-11 shows properties of waves 
at 3.^1 X 10«-3 Hz are similar to those of the 1.75 X 10”-4 Hz waves. 
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5.2) DIGITAL FILTERS (OPERATORS FILDES, FILOPT, FILTER AND VfTS) 

5.2.1) INTRODUCTION TO DIGITAL FILTER DESIGN USING REMEZ EXCHANGE 
Built into IDSP are operators that allow a user to design nonrecurslve Finite 
Inpulse Resfx>nse (FIR) digital filters and then apply these filters to the 
time series. Ihere is currently no capability to design Infinite Impulse 
Response (HR) filters in IDSP. 

Nonrecursive FIR filters have the following properties: 

1) They are implemented in a nonrecursive fashion (operator FILTER), 
i.e., there is no feedback in the equation, (thus the name "Finite"; 
HR filters must be implemented recursively) . 

2) Because of property (1) above, they are always stable. On the other 
hand the stability of IIR filters must always be checked. 

3) They can be made to be linear phase operations, i.e., with zero phase 
distortion, by making the number, N, of filter coefficients odd and 
centering about sample nunber (N-D/2. HR filters are not of linear 
phase . 

4) They are straightforward to design (operators FILOPT and FILDES). 

5) FIR filters are not as efficient as IIR filter, i.e., for given 
design specifications HR will require fewer multiplications. This 
efficiency can be important when large quantities of data are to be 
filtered and/or the filter has a large nunber of coefficients, or In 
the case vrtiere the filter is to be implemented in hardware. However, 
in IDSP these considerations are usually not important, so only FIR 
design is provided. 

The FILOPT and FILDES operators use a computer program developed by McClellan, 
Parks and' Rablner (see Rablner and Gold, 1975) that employs an algorithm 
called Remez exchange which is very efficient in producing an optimal (in the 
weighted Chebyshev sense) filter. This program has found wide application in 
the area of digital filtering and is considered state-of-the-art. The 
original program was obtained from the IEEE Library of Program*) for Digital 
Signal Processing (see Programs for Digital Signal Processing, 1979) and 
restructured as operators for use in IDSP. 


Before illustrating the use of these operators, we Introduce the coacept of 
normalized frequency: NORMALIZED FREQUENCY=( ACTUAL FREQUENCY) /(SAMPLE RATE), 
with units of cycles per sample. Thus the Nyquist frequency will always be 
0.5 in normalized Hz. As an example, if we sample a 20 Hz signal 100 
times/sec the normalized frequency is 0.2 Hz. 
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1 Tile important filter design specifications are shown in Fig. 5.2-1 (see 

*1 Schaff, 1979) where the abscissa is plotted as normalized frequency. 

•i 

I [0,f 1 Is the normalized PASSBAND 

[f ,f ] is the normalized transition width, delta _ 
p s r 

[f ,0.5] is the normalized STOPBAND 
s 

delta ^ is defined as the amplitude deviation in the PASSBAND 
delta^ is defined as the amplitude deviation in the STOPBAND 
PASSBAND RIPPLE is defined to be 201og^Q(1+delta^) in dB 
STOPBAND RIPPLE is defined to be -POlog^^Cdelta^) in dB 

..i Typical values for the PASSBAND RIPPLE would be in the range of 1 to 0.01 dB 

and for the STOPBAND RIPPLE 20 to 90 dB. 




5.2.2) DESIGN OF BANDPASS/BANDSTOP FILTERS 

The following example illustrates how to go about designing a filter and in 
turn applying it to the time series. Ori the Dynamics Explorer Spacecraft 
there is a magnetic field experiment that provides a vector measurement of the 
ambient magnetic field every 0.5 sec. The data from one of the components of 
the vector is shown in Fig. 5.2-2 as a function of time. As is evident there 
appears to be a substantial periodic signal riding on top of the actual signal 
of interest. It was desired to remove this contamination from the signal of 
interest. First it was necessary to do a spectrum analysis on the data to 
determine the frequency(s) of the unwanted signals present in the data. 

Fig 5.2-3 shows the results of using operators WINDOW, FFT, SPECT and GRAPH 
on this time series; note that at 0.0324 and 0.175 Hz there are peaks in the 
spectrin resulting from the unwanted signals. 

The first order of business is to design a FIR filter to remove the unwanted 
0.0324 Hz signal (0.0162 normalized Hz, e.g., 0.0324/2). 
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Referring to Fig. 5.2-1, once the designer has decided on the nunber of filter 
coefficients and the normalized transition width for a particular filter 
design the important parameters are 1 * delta ^ (the bandpass ripple) and 
deltag (the bandstop ripple) as a function of the weight ratio, delta^/delta^ 
s WTX(1)/WTX(2) (see Rabiner and (k>ld, 1975, PP 197,199). Delta^ is a measure 
of how successfully we will attenuate the unwanted part of the frequency 
spectrum and we want to make this value as small as possible, while at the 
same time minimizing the bandpass ripple, 1 + delta In order to aid in 
making the necessary tradeoff decisions we have developed an operator, FILOPT, 
which computes, using the Remez exchange algorithm, the values, in dB, for 1 * 
delta ^ and del tag over a range of delta^/deltag ratios. Typically operator 
FILOPT (see Section 6) would be used to compute 1 + delta ^ and delta^ over M 
decades of the ratio (from 0.1 to 1000). Because it is necessary to execute 
the Remez exchange algorithm over such a wide range of ratios it was decided 
to only do the computation for 9 values each decade and use operator INTERP to 
interpolate for the equidistant spacing required by IDSP. This dramatically 
decreases the amount of computation time required to execute FILOPT. The 
operator FILOPT produces an IDSP compatible data set consisting of elements of 
three time series for each weight ratio (see Fig. 5.2-4): series 1 s passband 
ripple = 201og( 1+delta^); series 2 s stopband ripple = 201og(deltag) ; and 
series 3 = the minimum value of delta^ in the stopband. It is then necessary 
for the designer to interpolate (operator INTERP) data set FILOPT so as to 
produce the data set, INTERP, that can then be graphically displayed using 
operator GRAPH. The designer can then use these graphs to make a decision 
about what weight ratio will give the best stopband attenuation for the 
minimum of passband ripple. 

Another approach, if there is some flexibility in the number of filter 
coefficients that can be used, and the designer wishes to prespecify the 
values for delta delta^ and the transition width is to estimate the nvnber 
of filter 'Coefficients from the following empirical formula which is valid for 
lowpass filters and has been Implemented as operator WTS (see Section 6). 
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NPILT « D(delta^.delta 2 )/(delta F) . f (delta ^.delta^) (delta F) -i- 1.0 

delta F = ABS(f.-f ) 

a p 2 

D(delta^.dalta2) » Ca^(log^.^delta^) a^log^^delta^ * a^]log^Qdelta2 

♦(a^(log^pdelta^)^ ♦ a^log^^delta^ ♦ a^] 
f(delta^.delta2) s b2(log^Qdelta^ - log^Qdelta2> 

a^a5.309X^0"^; a2«7.11*»X10”^; a^*-4.76lX10“^ a,s-2.66X10”^; a5=-5.9‘HX10”^ 
ag*-9.278X10"’; b^*11.01217; b2*0.512M 

In our case we used FILOPT to determine the best tradeoff for delta ^ and 

delta^. For this filter we decided to use 149 coefficients; fp s 0.0055 and 

f s 0.0155 for the first transition width and f = 0.017 and f = 0.027 for 
s s p 

the second transition width. We then use FILOPT to provide us with guidance 
on what is the best tradeoff for delta ^ and deltag. Referring to the output 
of FILOPT shown in Fig. 5.2-4 we conclude that a ratio of WTX(1)/WTX(2) = 2.5 
gives us small bandpass ripple and acceptable bandstop characteristics. We 
now have all of the design specifications for this first filter and can set up 
specifications to allow operator FILDES to be used to design the desired 
filter. 

To do this we generate a data set on disk that has the following format: 

First record: NFILT, JTYPE.NBANDS.LGRID. (415 format) 

Next record(s): EDGE (4F15.9 format) 

Next record(s): FX (4F15.9 format) 

Next record(s): WTX (4F15.9 format) 

These records can be placed into [usrid.IDSPlFOR058.DAT or a data set of 
another name, if the user does not want to enter design specifications 
interaotively. In general it is probably better to generate the data set for 
purposes of documenting the design rather than entering them interactively. 

NFILT— FILTER LENGTH; 

JTYPE— TYPE OF FILTER 

1 s MULTIPLE PASSBAND/STOPBAND FILTER, 

2 « DIFFERENTIATOR. 

3 « HILBERT TRANSFORM FILTER; 
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NBANDS-. NUMBER OF BANDS; 

LGRID— GRID DENSITY, WILL BE SET TO 16. 

EDGE(2»NBANDS)— BANDEDGE ARRAY, LOWER AND UPPER EDGES FOR EACH BAND WITH A 
MAXIMUM OF 10 BANDS. 

FX(NBANDS)— IDEALIZED PERFORMANCE ARRAY (OR DESIRED SLOPE It A 
DIFFERENTIATOR) FOR EACH BAND. 

WTX(NBANDS)— WEIGHT FUNCTION ARRAY IN EACH BAND. FOR A DIFFERENTIATOR, THE 
WEIGHT FUNCTION IS INVERSELY PROPORTIONAL TO F. 

In the case of this filter: 

NFILTr 1«9; filter coefficients 

JTYPE: 1; multiple PASSBAND/STOPBAND filter 

NBAND: 3; one STOPBAND and two PASSBANDS 

LGRID: 0; for purposes of plotting the transfer 

function-always set to 16 by the operator 

BANDEDGE ARRAY (in normalized freq): [0.0,0.0055]; 

[0.0155,0.0171; [0.027,0.5] Note that two records are 
required for the BANDEDGE ARRAY in this example because of 
the 4F15.9 format, and we have 6 nunbers to enter. 

FTX: 1.0, 0.0, 1.0; Idealized performance of the filter in each 
of the bands; in other words we wish no attenuation in the 
PASSBANDS and complete attenuation in the STOPBAND. 

WTX= 1. 0,2.5, 1.0; the ratio of delta^/deltag. 

Now that we have developed the design specifications for the filter, we can 
enter the specifications into a data set or enter them interactively; IDSP 
will do it either way. An example of the design specification data set is 
shown below. We then execute the operator FILDES, and it will produce a 
standard IDSP data set that contains the computed filter coefficients. 


149.1.3,0 

0.0,0.0055.0.0155.0.017 
0.027,0.5 
1 . 0 , 0 . 0 , 1.0 
1 . 0 , 2 . 5 , 1.0 
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Looking at Table 5«2-1, which is also produced by the execution of the 
operator FILDES, we see the coefficients and design specifications printed 
out. "LOWER BANDEDGE & UPPER BANDEDGE" gives us the bandedge array values; 
"DESIRED VALUE" gives the idealized performance; "WEIGHTING" gives the ratio 
of delta ^/delta 2 and "DEVIATION A DEVIATION in dB" give the amplitude 
deviation resulting from the design, 1.0 delta^ and delta^* Note Table 5.2-1 
only shows 63 coefficients because of symmetry. 

The next step is to look at the frequency response (transfer function) of the 
filter that was Just designed to ensure that it is what we want. To do that 
remember that the FIR filter is simply a "black box" through tdiich we are 
going to pass our time series, and what we want to unow is the frequency 
response of this "black box" for all frequencies from DC up to the Nyquist 
frequency (0.5 normalized frequency). The frequency response of a LINEAR 
SHIFT INVARIANT SYSTEM, in this case our FIR filter, can be obtained by taking 
the Fourier transform of the impulse response of the linear system. All 
linear systems produce an output by convolving the input signal with the 
impulse response. In our case the input signal is our time series and the 
Just calculated^ filter coefficients are the impulse response. If we take the 
Fourier transform of the filter coefficients we will obtain the frequency 
response of the filter. Note that it is only necessary to look at the 
amplitude part of the frequency response because we know by definition of an 
FIR filter that the phase is linear. 

After execution of operator WINDOW (using the rectangular window and zero 
padding options to Increase resolution of the plot) on data set FILDES (which 
contains the filter coefficients) , we execute the operator FFT with the 
magnitude/phase option (which stores the results in magnitude/phase form) on 
data set WINDOW. We then execute the operator GRAPH on data set FFTNP to 
produce Fig. 5.2-5 which is the amplitude response for the filter from DC up 
to 0.5 normalized Kz. Note that we used operator WINDOW and zero padded to 
better define the transfer function for purposes of plotting. Fig. 5.2-6 
shovo the original time series after being filtered with the FIR filter that 
we Just designed. 
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In a similar fashion we design the second filter, again using FILOPT to 
determine the optimal design, to remove the unwanted 0.175 Hz (or 0.0875 
normalized Hz) signal. The amplitude response of this second filter is shown 
in Fig 5.2-7. The design specification data set appears below: 

125.1.3.0 

0.0,0.05.0.075.0.1 
0.125.0.5 
1.0, 0.0. 1.0 
1.0,38.0,1.0 


Finally, we filter the time series a second time to remove the 0.175 Hz signal 
and the resulting time series is shown in Fig. 5.2-8, substantially free of 
both unwanted signals. 

It should be noted that for purposes of Illustration we have designed two 
separate filters; however, both unwanted frequencies could have been removed 
with a single filter that has two stopbands and three passbands for a total of 
5 bands. Operator FILOPT and FILDES will permit you to design a filter with 
up to a total of 10 bands and 512 coefficients. 

The design tradeoff information of this 5 band filter from FILOPT is shown in 
Fig. 5.2-9. The transfer function is shown in Fig. 5.2-10 where a 
delta^/delta^r 8.0 was chosen as optimal. Table 2.0-1 results from the 
execution of operator FILDES and presents the coefficients and other 
parameters of this 5 band filter. Fig. 5.2-11 shows the results of filtering 
the original time series from Fig. 5.2-2 with this filter. Fig. 5.2-12 shows 
the power spectrum of this filtered data processed with the 5-band filter; 
Fig. 5.2-13 is the power spectrin of the data set filtered with the two 3-band 
filter. In both cases the contaimination signals have been removed. 

Fig. 5.2-21 depicts the aeries of operations Involved in t^.e optimal design oi 
a multiple band filter using IDSP. 
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5.2.3) USE OP A PIPFERBNTIATING FILTER; 

The spectral analysis of a finite duration time series nay Include errors due 
to "leakage**. This error Is most likely to modify regions of the spectrum 
which contain little power by leaking power from more energetic frequencies, 
through the Influence of the side lobes of the data window (see Section 5.3 
for a detail discussion of windows) used In the analysis. One way to 

determine If a particular spectral feature is real or due to leakage Is to 
"pre-whlten” the series using a differentiating filter. In effect this 

multiplies the original spectrum by the "one" power of the frequency (F) which 
amplifies the high frequency relative to the low. For spectra which are 
dominated by low frequency power, the differentiator flattens, or "whitens" 

the spectrum and decreases the leakage effect on high frequency spectral 
features . 

As an example of designing a differentiating filter consider the tangential 

(T) and normal (H) components of a magnetl'* field vector plotted as a function 

of time In Fig 5.2**14 and Fig. S.2-15, respectively. Fig 5.2-16 shows the 

-5/3 

power spectrum of this vector. This spectrum decreases at approximately F 
for large P, so that spectral features In the higher F region may be 
Influenced by leakage from the lower F region. 

A differentiating filter was designed using the operator FIT-OES to prewhlten 

these data. This filter Is 126 coefficients In length with a cutoff frequency 

at 0.5 and a slope of 1. The design specification data s. t Is as follow: 
126,2,1,16 
0. 0,0.5 
1.0 
1.0 

The filter coefficients are shown In Table 5.2-2, and the .transfer function Is 
plotted In Fig. 5.2-17 The filtered versions of T and N are shown In Fig. 
5.2-18 and 5.2-19, respectively. The low frequency features are visibly 
absent, and the data appear to be much more like white noise. This Is borne 
out by the power spectrum shown In Fig. 5.2-20 which Is now flat at large F. 
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Fig. 5.2-1 

FIR filter design specifications after Schafft 1Q79* 
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Fig. 5.2-2 

Raw magnetic field time series from the DE spacecraft sampled at twice/sec, 
showing contamination riding on top of the desired signal. 



INDEPENDENT vfiPlfiBLE 

Fig. 5.2-3 

Power spectrum of the raw time series shown in Fig. 5.2-2 showing power at 
0.032'» and 0.175 Hz. 
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F1«. 5.2-4 relative weight 

The output of operator FILOPT showing 201og(1 ♦ delta.)* and 201og(delta2) 
plotced against the weight ratio* delta ^/deltag. From these plots w€ 
concluded that the best tradeoff was a weighi ratio of 2.5* which in turn waa 
Uv'ed to design the filter shown in Fig. -5.2-5. 
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Fig. 5.2-5 
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Transfer function of a notch filter designed to remove the unwanted 0.032M Hz 
signal from the original time series. 
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Fig. 5.2-6 

Time series after applying the filter shown in Fig. 5.2-5. 
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Fig. 5.2-7 

Transfer function of a notch filter designed to remove the unwanted 0.175 Kz 
signal from the time series. 
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Fig. 5.2-8 

Time series after application of both the filter shown in Fig. 5.2-5 end Fig. 
5.2-7. Note that the time series is substantially free of the unwanted 

signals. 
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The transfer function of the 5 band filter is shown. This transfer function 
is obtained by using operator FFT on data set WINDOW which contains the filter 
coefficients . 



Fig. 5.2-11 



INOEPuNDENT VRRIflBI.E 


Shows the results of using operator FILTER to filter the original time se.ies 
shown in Fig. 5.2-2 with the 5 band filter shown in Fig. 5,2-10. 
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Shows the power spectrun of the data after being filtered with the 5-band 
filter, note that the peaks in the power spectrum at 0.0324 Hz and 0.175 Hz 
have been removed. 
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Fig. 5.2-13 INDEPENDENT VRRIflBLE 

Shows the power speotrun of the data after being filtered with the two 3-band 
filters, note that the peaks in the power speotrum at 0.0324 Hz and 0.175 Hz 
have been removed. 
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Fig. 5.2-14 


Shows the tangential (T) component of a magnetic field vector plotted as a 
function of time. 
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Fig. 5.2-15 

Shows the normal (N) component of a magnetic field vector plotted as a 
function of time. 
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Fig. 5.2-16 

Shovfs the (,?wer spectrum of the magnetic field vector plotted in Fig. 5 
and 5.2-15. 
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Fig. 5.2-17 

Shows the transfer function of the differentiating flltrr. 
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Fig. 5.2-20 

The low frequency features are visibly absent and the data appears to be much 
more like white noise. This is borne out by uhe power spectrum shown in this 
plot which is flat at large F. 
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5.3) THE EFFECT OF WINDOWING. ZERO PADDING. AND OBSERVATION TIME ON THE 
ESTIMATION OF SPECTRA OF SINUSOIDS (OPERATORS SETUP. WINDOW. FFT, AND 
SPECT) 


There are two fundamental issues In estimating spectra (see operators WINDOW, 
FFT, AND SPECT) of sinusoids: 

Resolution- How close in frequency can multiple sinusoids be spaced and 
still be uniquely identified in the spectrum. 

Dynamic Range - How large can the differences in amplitudes of the 
sinusoids be and still be uniquely identified in the spectrum. 

The choice of the weighting function (window, see operator WINDOW) determines 
dynamic range and observation time (the length of the data record, see 
operator SETUP) determines resolution. 

The procedure is as follows: 

First we introduce the concept of normalized frequency which is defined as: 
NORMALIZED FREQUENCY= ( ACTUAL FREQUENCY) /(SAMPLE RATE), which has units of 
cycles per sample. So the Nyquist frequency will always be 0.5 in normalized 
Hz. 

a. ) Select weighting function (window) based on the desired dynamic range . 
See Table 5.3-1 below. The main lobe widths in Table 5.3-1 are given in 
normalized Hz. Reference 2 makes a strong case for using the 6 dB point 
as a performance indicator rather than the traditional 3 UB point. In 
other words if NslOO the main lobe width of the rectangular window is 
0.0121 normalized Hz. Note that all windows have a main lobe width 
proportional to 1/N normalized Hz. 

b. ) Select the number of data points, N, based on the desired resolution . 


TABLE 5.3-1 

WINDOW PEAK SIDELOBE MAINLOBE WIDTH 

(dB) (normalized Hz) 

(6 dB point) 


RECTANGULAR 

-13 

1.21/N 

HANNING 

-32 

2.00/N 

HAMMING 

-43 

1.81/N 


The following discussion explains how to go about this process : 
DYNAMIC RANGE (WINDOWING) 


An infinitely long pure sinusoid hdS as its Fourier transform a discrete line 
in the frequency domain. In practice, of course, we have a finite length data 
record of which we wish to take the Fourier transform. In effect we .see the 
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signal through a finite "window". If you do nothing but simply take N samples 
of a signal, you are windowing the signal ly uniform (rectangular) weights. 
In the frequency domain the spectrum of this data /ecord is the spectrum of 
this rectangular weighting function shifted by the sinusodial frequency. In 
the time domain you multiply the pure sinusoid by a rectangular window; in the 
frequency domain this becomes a convolution. 

For the "do nothing" (rectangular window) case the Fourier transform of this 
window Is plotted in Fig. 5.3-1 for a sinusoid with amplitude 50 and frequency 
20 Hz with one unit of noise added and a observation time of 0.5 sec. 

Note the 6 dB width of the main lobe is 1.21/N in normalized Hz. Therefore, 
if you do nothing but take N samples of a sinusoid and Fourier transform it 
you will get the function in Fig. 5.3-1 centered at the frequency of the 
sinusoid, which in this case is 20 Hz (0.2 normalized Hz'. 

As seen in Fig. 5.3-1 there is considerable energy in the side lobes. To 
reduce the size of the side lobes a different weighting function (window) 
whose side lobes are lower is used. However, the price paid for lower side 
lobes is an increase in the width of the main lobe. Using Table 5.3-1 and 
knowing the dynamic range of the signals of interest allows us to select the 
proper "window" . 

As an example suppose we are trying to resolve two sine waves with a ratio in 
amplitude of 1/20, further suppose in the frequency domain that the separation 
in frequency is such that the smaller (implitude sine wave falls at the peak of 
the first side lobe of the higher amplitude sine wave. From Table 5.3-1* if 
we use a. rectangular window, the side lobe will "mask" the smaller amplitude 
sine wave. Since a factor of 20 in amplitude is 201og.^(1/20)=-26 dD, the 
peak side lobe for the rectangular window at -13 dB will thus mask the smaller 
amplitude; however, the "banning" window is at -32 dB and thus the smaller 
signal should be visible. We have L.proved our dynamic range by a change of 
window! Note, however, that the 6 dB main lobe which was 1.21/N normalized Hz 
wide in the case of the rectangular window (see Table 5.3-1) has now widened 
to 2. 00/N normalized Hz with an attendant loss of resolvability . 

RESOLUTION 


Suppose there are two sinusoids present in the data record. In the frequency 
comain, after taking the Fourier ’.ransform (see operators FFT and SPECT), ire 
will have one or two peaks depending on how close in frequency they are 
relative to the main lobe of the particular window that was used. As the 
frequencies get closer together, the main lobes go from a bimodal to unimodal 
c urve . A commonly accepted criterion for resolution is that 1/2 of the main 

lobe width separation is needed to resolve two sinusoids. Using Table 5.3-1 
we can determine the main lobe width for the particular window used. 

EXAMPLE 


As 'n example of the complete process of selecting dynamic range and 
rpsulutlon, consider a signal sampled at 100 times per second, i.e., sample 
rate =100. Suppose two unknown frequencies are present in the data separated 
by at least 2 Hz, and we know that the amplitude of one frequency can be no 
greater than 50 times 'he amplitude of the other frequency. For this example 
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we have used a time series with two frequencies (20.5 Hz and 18 Hz or 0.205 
and C.18 normalized Hz) with amplitudes 50 and 1 units, respectively. We now 
need to answer the following questions: 

a. ) What is the dynamic range of the signals, thus what window should be 
employed? 

b. ) What is the desired resolution , thus determining the length of the data 
record, i.e, the observation time? 

The answer to question (a) is that the dynamic range is 201og.„(1/50)=-3i| dB; 
thus we need to use a Hamming window with a peaK side lobe is down at -43 dB. 
Fig. 5.3-2 shows the Fourier transform of such a time series which has been 
transformed using a "hanning” wirdow. Note the transform does not show the 
presence of the lower amplitude frequency. However, Fig. 5.3-3 is the same 
time series transformed using a "Hamming" window which clearly shows the lower 
amplitude frequency. Both Fig. 5.3-2 and 5.3-3 show the widening of the main 
lobes as the result of using the banning and Hamming window. 

The answer to question (b) is that the resolution necessary is 2.5 parts in 
100 or 0.025. Using the criterion that 1/2 of the main lobe width separation 
is needed to resolve 'wo sinusoids, from Table 5.3-1 the main lobe width c." 
the Hamming window is 1.81/N normalized Kz. So l/2( 1 .81/N)=0.025, and, 
solving for N, we get data points or a 0.36 sec observation time. Note in 
Fig. 5.3-4 with an observation time of 0.25 seconds we are not able to resolve 
the lower amplitude frequency even though a Hamming window was useo. Whereas 
in Fig 5.3-3 with a 0.5 sec observation time the lower signal is still seen. 

What this amounts to in practical terms is that if there are several features 
"close" together we must have a "longer" observation time or the convolving 
windows will not let us distingish them from each other. The "closer" the two 
features are to each other the "longer" the observation time req lired to 
resolve them. The resolution limit using Fourier ttohniquos and a **€ tangular 
window is 1/T^ where T^ represents the observation time. 

ZERO PADDING ;Zero padding is an option in Operator WINDOW' 

The purpose of zero padding is to better define, for purposes of appearance 
and accuracy, the plotting of the transform. Zero padding, as the name 
implies, is adding additional zeros to the time series before the Fourier 
transform is computed so that the transform wil be evaluated at a larger 
number of normalized frequencies than it would be if only the original N data 
points were used. 

When zero padding and using a rectangular window, one has to exercise caution 
because padding may introduce a discontinuity into the time series, i.e., the 
series can abruptly go to zero at the first padded point (see Fig. 5.3-5), 
this in turn will introduce additional unwanted spectral contributions 
(leakage). It should be noted that the Fourier transform implicitly assumes 
that the time domain function observed during T. is repeat<’d with period T. 
outside of the observation limits. Consequently when zero paddi.ig we have 
modified thr original time domain function and consqiiently would expect this 
modified function to have a different transform and no longer approximate the 
original continuous Fourier transform. 
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In Fig. 5.3-6 we show the transform of a 10 Hz time series, using a 
rectangular window, without zero padding and, i.. Fig. 5.3-7, the transform 
after padding; notice the much higher leakage levels introduced by the 
discontinuity due to the padding. Recall that windows are applied to reduce 
the order of the discontinuity at the boundary of the periodic extension, so 
in the case of a time series windowed with other than a rectangular window, 
zero padding will not introduce a discontinuity because the window tapers the 
time series to zero prior to padding. This tapering can be seen in Fig. 5.3-8 
where we show the time series of a 10 Hz signal after application of a Hamming 
window. Fig. 5.3-9 and 5.3-10 show the transform of a 10 Hz time series after 
application of a Hamming window without zero padding and with zero padding, 
respectively. Note that there is only a small increase in the leakage 
introduced by the padding. 

Note that zero padding does not improve resolution since only an increase in 
observation time (T^) will do that. But zero padding does result in 
additional Fourier coefficients being computed which better defines the 
sampled continuous Fourier transform. 
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Fig 5.3-1 

Transform of a sinusoid of amplitude 50 units, and frequency 20 Hz with one 
unit of random noise added using a rectangular ("do nothing") window. 
Observation time of 0.5 sec. 



‘ Transform of a time series containing two frequencies (20,5 and 18 Hz) with 

amplitudes 50 and 1 units, respectively, and one unit of random noise using a 
banning window. Observation time of 0,5 sec. Demonstrates the inability of 
the banning window to resolve the two frequencies present. Noce the width of 
the main lobe has increased relative to the main lobe width of the rectangular 
window. 
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INDEPENDENT VRRlfiBlE 

Fig. 5.3-3 • 

Transform of a time series containing two frequencies (20.5 and 18 Hz) with 
amplitudes 50 and 1 units, respectively, and one unit of random noise using a 
Hamming window. Observation time of 0.5 sec. Demonstrates the ability of the 
Hanning window to resolve the two discrete frequencies. 



/ independent vofiinniF 

j Fig. S.S-** 

I Transform of a time series containing two frequencies (20.5 and 18 Hz) with 

amplitudes 50 and 1 units respectively and one unit of random noise using a 
Hamming window. Observation time of 0.25 sec. Demonstrates the inability, 
even using the Hamming window, to resolve the two frequencies present due to 
inadequate observation tine. 
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Plot of a 10 Hz time series showing the discontinuity introduced by zero 
padding. The Fourier transform implicitly assumes that the time domain 
function observed during Tq is repeated with period Tq outside of these 
limits. Consquently zero padding will modify the original function and result 
in a different transform. 





Fig. 5.3-6 

Transform of a 10 Hz time series containing 1 unit of noise with no zero 
padding using a rectangular window. 
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Fig. 5.3-7 

Tr«nsforB of a 10 Hz time series containing 1 unit of noise zero padded to 200 
points using a rectangular window. Note the increase in leakage due to zero 
padding. 



Pig. 5.3-8 


Plot of a 10 Hz tine series after being windowed by a Hanning window, note how 
the ends are tapered to zero. 




Fig. 5.3-9 

Transform of a 10 Hz tiaw series containing 1 unit of noise using a Hanning 
window and no zero padding. 



Fig. 5.3-10 

Transform of a 10 Hz tlM series containing 1 unit of noise using a Hanning 
window zero padded out to 200 points, note that the leakage level is about the 
same as shown in Fig 5.3-9* 
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6.0 IDSP OPERATOR DESCRIPTIONS 


Usage of AVER Command 


The user specifies 'AVER dsname n* where n must be the number of points which 
he wants averaged into each box. Tne last box may not be full depending on how 
n divides the number of available points. If the user specifies optional 
parameter ZERO, the first box in the average contains only one point. (This 
option can be used when one wants to avoid having the zero frequency averaged 
in with anything.) The following exa.iple utilizes this option, with the rest 
of the points of dataset SPECTOFF averaged into 8 points per average: 


AVER SPECTOFF 8 ZERO 



Usage of CMDHIS Command 


This command causes a printout of user responses to the system query ’'ENTER 
COMMAND" that were entered during the current analysis session. It has no 
operands. The file that contains these commands is FOR022.DAT and can be 
printed after the session ends. 
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Usage of CONCAT Comnand 


If two datasets have the same nunber of points (NPTS), if both are of the same 
type (real or complex) « and if the sum of their series is less than or equal 
to eight, they can be concatenated into one dataset of NPTS points and (NSER1 
+ NSER2) series. The form of the command is; 

CONCAT datal dataZ 
For inverse operator, see SUBSER. 
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Usage of the COPY Comnand 


The IDSP system destroys datasets in its default working directory 
[usrid.IDSP] by overwriting the old dataset when a given operation is 
repeated. After any operation, any dataset can be copied into a file of some 
other name. The first example copies the current version of WINDOW into WIND6 
in the working directory; the second copies it into a different directory file 
[usr id. GAMMA ]WIND7. DAT. 

COPY WINDOW WIND6 

COPY WINDOW [ usr id. GAMMA 1WIND7 
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Usage of the DCL Conn and 


This operator allows the user to temporarily exit the IDSP process and execute 
VAX Connand Language Coimands via a subprocess. He will receive the standard 
prompt. When finished the user should type: "STOP/IDsO" or "LO" This will 
return control to IDSP by terminating the subprocess. NOTE: Typing 
produces a message indicating the subprocess was logged off. It is no '• ■ 
for alarm. The user is now back in IDSP. 

WARNING: Do not do something foolish in DCL like Control Y, or a delete, u. a 
purge. When in DCL your LOGIN constraints, such as "decide" for delete o. 
KEEP=2 for purge are not recognized. 
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Usage of DSTAT Command 


To obtain statistics of a dataset, type "DSTAT dataset- name". Any data 
flagged as BAD (see SETUP) will not be included in the statistical 
calculation; the percent of bad data per series will be calculated. Statistics 
calculated for each series: 


0 

XM 

Mean 

rx^/n 

1 

XMABS 

Mean of absolute values 

zlX^I/n 

2 

XMSQ 

Mean of values squared 

£X^^/n 

3 

XRNS 

Root Mean Square 

SORT (XMSQ) 

4 

VAR 

Ifribiased variance (division by n-1) 

£(Xj-XM)^/(n-1) 

5 

NBAD 

Number of bad data per series 


6 

PER 

Percent bad data per aeries 
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Usage of the EDHTST Comnan<i 


Ihe EDHIST operator allows the user to edit the history elements in the header 
of j specified dataset., The form of the command is; EDHIST dsname where 
"dsname" is the name of the input dataset. 

This operator has several editing commands listed below: 

D: Delete the current line. 

R: Replace the current line. 

I; Insert a line before the current line. 

<cr>; Save the current line. 

E: End of editing session. 
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Usage of EDIl ommand 


The interactive user at a HP-?648A Graphics Terminal can dynamically flag data 
points as BAD (set equal to -99«9) in any real dataset. The EDIT operator 
causes a display of up to 500 points to the screen and places the user into 
interactive usage of interactive graphics comnands. (To see a list of these 
commands at this stage, type "HELP".) The only commands that have any effect 
on the dataset stored by EDIT are those that perform deletions — the EDIT 
operator will flag them as BAD as soon as the user types "EXIT" to return to 
the EDIT operator. If the series had more than 500 points, the next group of 
points will then be displayed. As soon as editing is finished on one series, 
EDIT moves to the next one. The form of the command is: 

EDIT dataset SEH=J ALL 

where "SERsj" is optional if the user wants only to edit the jth series. "ALL" 
is optional and specifies that if the ith point is deleted in any series, then 
it is deleted in all the series. 

Hore extensive documentation on the HP interactive graphics commands can be 
found in Carleton, T. P.. LEP VAX Users Guide. NASA TM 84931# October 1982. 
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Usage of the EIG Comnand 


^7 


■ « 


I 

I 


This operator takes as Input spectral matrices and does rotation to an 
eigenvalue coordinate system (see Section 5.1 for details). Statistics for 
some of the computations are extracted from a dataset created previously by 
operator DSTAT, operating on a time series data set such as the output of 
INTERP or FILTER, i.e., a pre-FFT data set. The form of the command is: 


EIG SPECTD SPECTOFF DSTAT MF IKON 


where "SPECTD" is the dataset containing the diagonal terms of the spectral 
matrix; "SPECTOFF" contains the corresponding off^Jlagonal terms; and "DSTAT" 
contains the dataset of series statistics. "MF" is optional and is used only 
if the user wants the analysis to assume that the original series have been 
placed into a Mean-Field coordinate system (see operator MNFLD). The diagonals 
of the rotated system are stored in EIGD, the off-diagonals are stored in 
EIGOFF, and various scaler parameters are stored in EIGPARM, and the 
eigenvectors and eigenvalues associated with the minlmun, intermediate and 
maxiraim directions are stored in EIGVEC and EIGVEC1. 

KCON is an optional parameter of EIG which can be used to select a sign 
convention for the wave normal vector li determined in the polarization 
analysis (but with arbitrary sign). KCON can have the values 1, -1 and 0, 
with KCON = 1 the default state. In this case, the component k. is forced to 
be positive by inversion of 1c (multiplying all components by -T) if found to 
be negative. For coordinate systems centered at the sun, that will force 1c to 
have a component outward (along the radius vector direction) from the sun. 
For data coordinate systems defined with X positive toward the sun, the user 
should choose KCON = -1; this forces k^ to be in the negative X direction by 
Inversion of 1c, again insuring that 1c have a component outward from the sun. 
If KCON = 0 is chosen, a different convention, based on the direction of 
average is applied. In this case the quantity BDOTK will be tested; if 
BDOTK < n, then ft is inverted and BDOTK reccxnputed. This forces ft to have a 
component in the +1 direction. 

Data set EIGPARM contains scaler quantities DEGPOL, COH, ELIP, THETA, TRTOT, 
SNR, BDOTK, and SIGMA. EIG also creates data set EIGVEC which contains in 
series 1-3, the first three components of the eigenvector corresponding to the 
direction of minimum variance for that mode. Series 4-6 are the eigenvalues 
(minimun, intermediate and maximim, respectively). Also created by operator 
EIG is data set EIGVEC 1 which contains the three components of the eigenvector 
corresponding to the intermediate direction of variance (series 1-3) and the 
three eigenvectors corresponding to the maximum direction of variation (series 
4-6). 

WARNING: If the "MF" option has been specified, the normalized magnetic 
hellcity (SIGMA) in EIGPARM will not be correct. 

For a schematic of operator EIG and its relationships to operators DSTAT and 
SPECT and the associated data sets see Fig. 5.1-8. 
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Usage of FFT Conmand 


The IHSL routine FFTRC is used to compute the Fourier Transform, X, of every 
series. A, in the dataset according to the following formula: 

X(K-4.1) = SUM FROM J = 0 TO N-1 

OF A(J+1)«CEXP((0.0,(2.0«PI«J«K)/N» FOR K=0, 1, . . . ,N/2 AND PI=3.1415... 

Only N/2+1 points are stored because of synsnetry. (If the number of points in 
each input series was odd, it is forced to be even by deleting the final point 
of each series before taking the transform.) The default normalization factor 
is 1. The user can choose his own normalization by specifying the optional 
.arameter: NORM=x where x is the value he desires. The first of the following 
examples has default noraalization; the second has normalization .000005. 

FFT INTERP 

FFT WINDOW N0RM=5.E-6 

WARNING: If N/2»2 has a large prime factor, the "FFP* will be VERY slow. 

If the user wishes to obtain the magnitude and phase of the FFT, he can 
specify the optional parameter "MP” which will also produce a complex dataset 
FFTMP where the magnitude is stored as the real part, and the phase in degrees 
is stored as the imaginary part. 

For inverse transform, see FFTIN. 

LIMITATION 


Because of the size of working arrays required, operations involving the FFT 
are limited to a much smaller size than other operations (such as SETUP). For 
series involving more than 10,000 points, the user should check docunentation 
under "DESIGN CONSIDERATIONS" and also be aware that the time for such long 
runs may be excessive. 
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Usage of FFTIN Conmand 


The IMSL routine FFTCC is used to compute the inverse Fourier Transform, X, of 
every series, A, in the dataset, each containing NPTS. First, a working array 
is constructed containing the complete transform, using the conjugate symmetry 
of the data: 


N = (NPTS-1)»2 

B(N+2-I) = CONJG(A(D), FOR 1=2, ...,N/2 
Then X is computed according the the following formula: 

X(K+1 ) = SUM FROM J = 0 TO N-1 

OF B(J+1 )»CEXP((C.0,(-2.0»PI»J*K)/N)) FOR K=0, 1, . . . ,N-1 AND PI=3.1415... 

The default normalization factor is (1/N). The user can choose his own 
normalization by specifying the optional parameter: NORM=x where x is the 
value he desires. The first of the following examples has default 
normalization; the second has normalization .000005. 

FFTIN SPECTOFF 

FFTIN SPECTD NORM=5.E-6 

For the forward transform, see FFT. 

LIMITATION— see FFT 
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Usage of FILOES 


For data having sample spacing T, the highest observable frequency (called the 
Nyqulst frequency) Is FN = 1/(2*T). To formulate filter design parameters, use 
a normalized frequency scale S = [0.0,0.53 where 0.5 corresponds to the 
Nyqulst frequency (FM). Thus, given any frequency F that the user wishes to 
input as a filter design specification, he first scales It by the formula: 

FS s 0.5 • (F / FN) or equivalently 

FS s F • T 

Hence, If the user wishes to enter a ft*equency parameter corresponding to a 
frequency of 20 (as in cos(2*PI*20*t)), and If his data spacing is T : 0.01, 
then FN s 50, and FS s 0.2. 

To execute the command, type a sequence of the form "FILDES r DSNsname** where 
"r** Is a required parameter specifying the source device (T s terminal, F = 
disk file) for the filter design specifications. When reading from a file, the 
user can specify the optional parameter DSN to equal the name of the 
particular file he wishes to access. When this parameter is omitted, the 
latest version of F0R058.DAT is always used. The following example executes 
the design according to the fourth version of the design spec file: 

FILDES F DSN sL0WPASS1.DAT; 4 

The filter coefficients will be stored in a dataset called FILDES. Ihe 
interactive user who desires to see the effect of several different filters 
should copy FILDES into some other dataset name (see COPY) before executing 
FILDES again. Then he can (as an option to the FILTER command) choose which 
filter he wants applied to his data. 

Filter design by Remez exchange technique. Maximum filter length is 512. Input 
data records as shown below: 

First record: NFILT, JTYPE,NBANDS,LGRID. (415 format) 

Next record(s):EDGE (4F15.9 format) 

Next record (s): FX (4F15.9 format) 

Next record(s): WTX (4F15.9 format) 

These records can be placed into [usrld.IDSP3FOR058.DAT if the user does not 
want to enter design specifications interactively. 

NFILT~ FILTER LENGTH 

JTYPE— TYPE OF FILTER 

1 s MULTIPLE PASSBAND/STOPBAND FILTER 

2 8 DIFFERENTIATOR 

3 * HILBERT TRANSFORM FILTER 

NBANDS— NUMBER OF BANDS (LIMIT OF 10) 
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LGRID— GRID DENSITY, WILL BE SET TO 16 UNLESS SPECIFIED OTHERWISE BY A 
POSITIVE CONSTANT. 

EDGE(2*NBANDS)— BANDEDGE ARRAY, LOWER AND UPPER EDGES FOR EACH BAND WITH A 
MAXIMUM OF 10 BANDS. 

FX(NBANDS)-- DESIRED FUNCTION ARRAY (OR DESIRED SLOPE IF A DIFFERENTIATOR) FOR 
EACH BAND. 

WTX(NBANDS)-- WEIGHT FUNCTION ARRAY IN EACH BAND. FOR A DIFFERENTIATOR, THE 
WEIGHT FUNCTION IS INVERSELY PROPORTIONAL TO F. 

SAMPLE INPUT DATA SETUP; 

32 , 1 , 3.0 

0.0,0.1,0.2,0.35 

0.425,0.5 (IN 4F15.9 FORMAT) 

O.C, 1.0,0.0 

10.0,1.0,10.0 


THIS DATA SPECIFIES A LENGTH 32 BANDPASS FILTER WITH STOPBANDS 0 TO 0.1 AND 
0.425 TO 0.5« AND PASSBAND FROM 0.2 TO 0.35 WITH WEIGHTING OF 10 IN THE 
STOPBANDS AND 1 IN THE PASSBAND. THE GRID DENSITY DEFAULTS TO 16. 

THE FOLLOWING INPUT DATA SPECIFIES A LENGTH 32 FULLBAND DIFFERENTIATOR WITH 
SLOPE 1 AND WEIGHTING OF 1/F. THE GRID DENSITY WILL BE SET TO 20. 


32,2,1,20 

0,0.5 

1.0 

1.0 


If the user wishes to peruse the printout describing the filter he can route 
the output* to disk via the FLOP conmand before executing the FILDES comnand. 

LIMITATIONS ; Avoid designing a differentiator with an ODD number of 
coefficients. At times the design algorithm complains about possible 
convergence problems. When designing a filter, ALWAYS MAKE A PLOT OF THE 
TRANSFER FUNCTION to ensure that the filter behaves as you desire. Section 
5.1.5. and particularly Section 5.2.2, discuss a technique for the systematic 
determination of the stopband ripple and passband ripple using operator 
FILOPT. Number of bands is limited to 10 for bandpass/bandstop filter. The 
maximum mmber of coefficients is 512 for any type of filter. 
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ORIGINAL PAG:- 
OF POOR QUALITY 

To obtain a plot of the transfer function of the filter, window the filter 
coefficients with a rectangular window and pad with enough zeroes to isolate 
the function. Specify: 

WINDOW FILDES 0 PADr512 

FFT WINDOW MP 

This creates a complex dataset FFTMP where the magnitude is stored as the real 
part, the phase (in degrees) as the Imaginary part. FFIMP can then be plotted 
using the GRAPH operator. 

More information available in Chapter 5 of IEEE PROGRAMS FOR DIGITAL SIGNAL 
PROCESSING (NASA/GSFC Library: QA 76.5.P77). See also Rabiner and Gold, THEORY 
AND APPLICATIONS OF DIGITAL SIGNAL PROCESSING, Chapter 3. 
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Usage of FILOPT 

This operator assists the user in determining the optimum ratio of delta ^ to 
delta 2 (See Fig. 5.2-1 for definitions of delta ^ and delta 2 ) to use in 
designing a filter with the Remez exchange algorithm. 

To execute the command, type a sequence of the form "FILOPT r 0SN=name" where 
"r" is a required parameter specifying the source device (Tsterminal, F=disk 
file) for the filter design specifications. When reading from a file, the user 
can specify the optional parameter DSN to equal the name of the particular 
file he wishes to access. When this parameter is omitted, the latest version 
of FOR059.DAT is always used. The following example executes the algorithm 
according to the fifth version of the design specification file: 

FILOPT F DSN sLOWPASS.DAT; 5 

Output is stored in a data set called FILOPT. THis file contains the following 
three series of information: 

1 - pass band ripple in db - defined as 20*log(1+delta1 ) 

2 - stop band ripple in db - defined as 20 *log(delta 2 ) 

3 - lowest limit of the stop band 

The values are calculated nine tiroes per decade for a maximum of five decades. 
The resolution of the output file is equal to the value of the smallest 
decade. All intermediate points are written as fill data (-99*9) to allow for 
subsequent interpolation in order to meet the equidistance spacing requirement 
of lOSP compatible data sets. 

Each filter is designed by the Remez exchange techique. Maximum filter length 
is 512. Input data records as shown below: 

First record: NFILT,NBANDS,LGRID 

Next reoord(s): EDGE (4F15.9 format) 

Next record (a): FX (4f15.9 format) 

Last record: IDEC,NDEC 

These records can be placed into [usrld.ldsp3FOR059.DAT if the user does not 
want to enter design specifications interactively. 

NFILT - filter length 
NBANDS - number of bands 

LGRID - grid density, will be set to 16 unless specified otherwise by 
a positive constant less than 16 

EDGE (2*nbands) - bandedge array, lower and upper edges for each band 
with a maximum of 10 bands 

FX(nbands) - desired function array for each band 

IDEC - starting decade for calculation (power of 10) 

NDEC - number of decades (maximum of five) 
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SAMPLE INPUT DATA SETUP: 
32,3,0 

0.0,0. 1,0.2, 0. 35 
0.425,0.5 
0.0, 1.0, 0.0 
-1,4 


This data specifies a length 32 bandpass filter with stopbands 0.0 to 0.1 and 
0.425 to 0.5, and passband from 0.2 to 0.35. The weighting ratio starts with 
the value 0.1 and ends with the value 1000.0 for the stopbands and a value 1.0 
in the passband. The grid density defaults ‘ j 16.0. 

LINITIONS: There are some problems in the filter design program that result in 
a lack of convergence message. The maximum number of coefficients is 512. 

To obtain a plot of the FILOPT dataset use operator INTERP to interpolate the 
fill data points and plot the values using the GRAPH operator. Specify: 

INTERP FILOPT 
GRAPH INTERP 


Warning: FILOPT operator executes the Remez exchange ^Igo^rlthm up to 50 times 
and creates a dataset containing three series of up to 100,000 points in 
length. This procedure takes a large amount of time and space to accomplish. 

Also, the GRAPH operate: uses a considerable amount of time to plot these 

values. For designs using four or five decades the run should be submitted as 

a batch Job. . |\ 

, V I V USER MAKES 

I USER X I 1 /" \ I 1 ^ DECISION ON 

SPECIFIES I OPR FILOPT! INTERP ^ akl am RATIO TO USE 

IK AF AND # ’ IFILOFTP ^INTERP W 'NTERP -W GRAPHjpr^^ ^^AND PREPARES 

I OF WEIGHTS I X. / I V-X J '' , X SPEC. FILE 


USER 
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FUNCTION 


ITERATE 
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Usage cf FILTER Command 


Subroutine FILTER extracts the filter coefficients from a dataset called 
FILDES unless the user specifies that they exist elsewhere by including the 
optional parameter ”DSNrdataset” • The dataset of filter coefficients must be 
created by FILDES before FILTER is executed. Filtering makes contiguous 
intervals noncontiguous. To obtain contiguous filtered intervals, use SETUP to 
create one large interval, filter it, and then use the SUBSET operator to 
produce the set of contiguous intervals desired. To filter the dataset WINDOW 
with the filter coefficients currently stored in FILDES, the user would 
specify: 


FILTER WINDOW 

If he wished to use filter coefficients stored in N91HP04 the user would 
specify : 


FILTER WINDOW DSN=N91HP0^ 

To avoid accidental aliasing problems, decimation can be accomplished only as 
an option of the FILTER command. The user should enter "FILTER dataset DECsn" 
where ”n" specifies that only the 0th, (n)th, (2n)th,. . . points will be kept 
after filtering. If this option is omitted, every point is kept (no 
decimation). WARhING: The decimation rate ”n" must be chosen so that l/(2*n*T) 
is greater than the lower edge of the lowpass filter stopband where T is the 
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Usage of FLOP Comnand 


Output for a given operation is routed according to the state of a system 
variable lOUT. The interactive user by default receives output at the 
terminal. If he wishes to have hard copy of some operation, he simply types 
the comnand "FLOP". From that point, printout goes to [usrid.IDSP3lDSP0UT.DAT 
which the user can access according to his needs after the session ends. 
Output to the terminal can again be restored by typing "FLOP". Interactive 
user prompts are unaffected by this command. 
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Usage of GRAFCK Coitnand 


This coRimand allows the user to check on the status of the subprocess that he 
initiated via the GRAPH command. The form of the command is: 

GRAFCK n 

where "n" is optional. If "n" is omitted, only one status check is made. 
Otherwise, a status check will be make every 15 seconds until either the 
subprocess has finished or "n" checks have been made. If the subprocess 
failed, a listing of the log file will be printed. Note that under some 
conditions GRAFCK will claim that the graph terminated normally even when the 
subprocess has actually failed. A look at the file (GRAF.LOG) will describe 
any problems in detail. 
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Usage of GRAPH Command 


Plot all the series of any dataset on the VERSATEC plotter in the specified 
type plot (0 = linear( x)-linear( y) , 1 = linear( x)-log( y) , 2 = 

log(x)-linear(y) , 3 = log(x)-log(y) ) . 

Specify line-type desired by (-1 = points plotted as symbols only, 0 = points 
connected by lines only, 1 = points plotted as symbols connected by lines). 

If only the j-th series is desired, it alone can be plotted by specifying the 
optional parameter SFR=J. 

Both real and imaginary parts of complex datasets are plotted on separate plot 
frames. The user can specify one or the other by the optional parameter "REAL" 
or "IMAG". 

The length of the x-axis can be set by specifying "XLEN=x" where x is in 
interval [1,100] (units of inches). The range of values desired in the plot 
can be specified via the option SIZE. If this option is chosen, the next 
record in the command stream must specify XL, XH, YL, YK in 4G10. 0. 

To obtain nice labeling in linear plots, specify limits so that the desired 
range Is neatly divisible by five. Logarithmic plots are forced to begin and 
end on decades. Points that fall outside the plot window are plotted at the 
boundary. The Independent variate.' is by default plotted as indices where i s 
0, 1, 2, ... , NPTS - 1. 

If the user wants it in terms of the actual delta spacing of the dataset, 
these indices will be multiplied by delta if the user specifies option 
"DELTA" . 

If plots are specified to be logarithmic in the first variable, the first 
point of the plot is automatically omitted. 

In the first example, a log-linear, symbols-only plot is given for every 
series of dataset WINDOW; in the second, only the third series of dataset 
INTERP is given a linear-linear, symbol-lins plot; in the third, only the real 
part of dataset FFT is given a linear-log, symbol plot with actual data 
spacing in the Independent variable and range limits for X and then for Y: 

GRAPH WINDOW 2 -1 

GRAPH INTERP 0 1 SERs3 

GRAPH FFT 1 -1 REAL DELTA SIZE 
0.,1.25, 10.,1.E5 
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If the user wishes to stack several graph coiAuands in rapid succession and 
then plot them all at once, he should specify the option "HOLD'* on all of them 
exceiit the last graph command executed. (WARNING: When this option is used a 
dataset name should usually appear only ONCE in the group being held for one 
graphing job. Remember that only the last version of the repeated dataset will 
exist at graphing execution. Depending on the command sequence this final 
version may have parameters incompatible with an earlier version, using the 
job to fail.) 

Special effects can be accomplished by combining other operators with the 
GRAPH command. For example, suppose the user wanted a logaritr.iiic plot of only 
the negative values of a dataset SPECTOFF. He could enter: 

NORM SPECTOFF -1 

GRAPH NORM SIZE 

where in the size parameter he would specify how closely he wished the plot to 
approach zero. The results that are "off scale" would correspond to 
non-negative values in the original dataset. 

The actual plotting algorithm accomplishes a VERSATEC plot via a subprocess 
which IDSP creates for the user. The first five characters of this process 
name will be the USRIO. To check the status of this subprocess use command 
GRAFCK. 

A log file of each graphing job is written to GRAFJ6.LOG so that the user can 
examine it in case his results were different from his expeccations. The user 
can easily delete all such files by the DCL command: "DEL GRAFJB.LOG.*". 

If the user has an HP2648A terminal, he can have graphics output sent to the 
screen by entering the option "DEV=HP". Up to 500 points of a given series 
will be displayed. Then the user can utilize the commands inside the 
HP-graphlcs package (for help, type "HELP" after points are displayed). When 
he has seen enough of those points, the user types "EXIT" to return him to the 
next segment of graphic display. He will then see his current graphics display 
plus an over-laid alphanumeric display describing the graphics. Graphic 
display can be turned on and off via key: "SHIFT GDSP"; alphanumeric via: 
"SHIFT ADSP". 

When the option DEV is usedj the following options are ignored: XLEN, HOLD. 
Also even though the linn-type parameter is required, the present 
configuration displays points only. (Inside the HP-package, these esn be 
connected by a line by entering "DATA LINE".) 
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Usage of the INTERP Cotimand 


Interpolation is linear , 
interpolated values. The 


If any flag values are found, they are replaced by 
ex2nple below results in linear interpolation: 

INTERP SETUP 
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Usage of MEM command 


Maximun Entropy Method of computing Povfer ^ectra (Anderson’s derivation^ 
Geophysics, Vol 39t Feb., 1974) 

This operator handles both real and complex datasets as input. The power 
spectrum is output as a real dataset (MEM. DAT). The subroutine COEFF 

calculates the filter coefficients A and the subroutine CPSPEC calculates the 
power spectrum at frequencies from "kero to f^ (Nyquist frequency) with the 
desired frequency spacing (FRQSPC). The default normalization factor (NORM) 
is 1. The default nunber of filter coefficients (NCOEF) is equal to 2N/ln2N 
(reference: Berryman, Geophysics Vol 43, Dec. 1978). The default nunber of 
frequency points (NFP) is equal to N/2 + 1, where N is the nunber of points in 
the original time series. 

'he optional paremeters may be input as follows: 

NORM sX; where X is the desired normalization factor (real) 

NCOEF sY; where Y is the nunber of filter coefficients (integer) 

NFD sZ; where Z is the nunber of frequency points (integer) 


EXAMPLE: 

MEM SETUP 

uses defaults values 

MEM SETUP NCOEF s50 NFP=75 
NCOEF = 50, NFP = 75, default normalization, NORM = 1 

MEM SETUP NORMs5.E-6 

uses given normalization factor and gets default values for NCOEF and NFP 


WARNING: The MEM operator is recommended for datasets with less than 1000 
points, increasing the nunber of points, NCOEF or NFP will cause MEM to be 
slower. 

LIMITATIONS: Because of the size of the working arrays required, operations 
involving MEM are limited to a much smaller size than some other operators. 
Therefore, datasets must contain less than 15,000 complex or real 
points/ series for up to 8 series. 



Usage of the MNFLD Coranand 


The operator MNFLD performs a transformation from the coordinate system 
of an input vector data set to mean field coordinates. It initially computes 
the mean field vector from the input data. In the new coordinate system, ZMF 
is taken in the direction of the mean field. XMF is chosen to lie in one of 
the original coordinate planes, determined by the Intended use of the 
transformed data. The desired plane is selected in response to a prompt: 1 = 

XY, 2 = YZ, 3 = XZ, which appears after the operator is commanded: MNFLD 

SETUP, MNFLD INTERP, etc. For whatever choice is given for XMF, YMF is 
properly oriented to complete the right-handed coordinate set. The 
transformation is only valid for matrix size s 3. 

If spectral analysis is to be performed it a later stage, note that 
transformation to mean field coordinates permits specification of the MF 
option on the subsequent call to the operator EIG, for example, 

EIG SPECTD SPFCTOFF DSTAT MF. 

This sets a flag MFIELD = 1 internal to EIG and its supporting subroutines, 
and enables application of a special eigenvector coordinate system definition 
appropriate for mean field data (see operator EIG documen«.atlon) . 

If the magnitude of the projection of the mean field vector on the plane 
chosen for the new X-axis is less than a specified lower limit (currently 
coded into subroutine MFCOMP, which is called by MNFLD, as the quantity BLO s 
0.01 times the total mean field magnitude), then the rotation matrix for the 
transformation is specified rather than computed (see matrix definitions 
following statement nunbers 61 and 81 in a listing of subroutine MFCOMP for 
the cases where XZ and YZ planes are selected, respectively). If the input 
data has a mean vector already in the ZMF direction (with an XY plane 
projection less than or equal to BLO), then no transformation is carried out, 
and the message ”N0 ROTATI'.J - DATA ALREADY IN MF COORD** appears. A listing 
of the mean field vector, and the transformation matrix used will appear on 
the terminal screen (or the above "NO ROTATION" message) unless a FLOP is 
commanded, which routes the output listings to IDSPOUT.DAT for later hardcopy 
display availability. 
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Usage of NORM Connand 


The NORM operator multiplies every point in the dataset by the specified 
factor. In the following example, dataset FFT is normalized by 0. 5E-4. 


NORM FFT .00005 



* 
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Usage of the RECPOL Cornnand 

The operator RECPOL will convert from rectangular to polar or polar to 
rectangular coordinates depending on parameter ITYPE. 

If ITYPEsI; rectangular to polar, if ITYPE=0; polar to rectangular. 

The data set created as output will be RECPOL.DAT. The location of x,y,z in 
the input data set will need to be supplied if converting from rectangular to 
polar, while the location of r, theta, phi in the input data set will need to 
be supplied if converting from polar to rectangular. If polar to rectangular 
the order of output data set will always be IXsx, IY=y, IZsz. If rectangular 
to polar the order of output data set will always be IX=r, lYstheta, IZsphi. 

EXAMPLE OF INVOKING RECPOL: 

RECPOL datasetname 


The user will then be prompted to enter whether he/she wants to go from polar 
to rectangular or rectangular to polar. ( ENTER A 0 IF WANT TO GO POLAR TO 
RECTANGUUR, ENTER A 1 IF WANT TO GO RECTANGUUR TO POUF. ) on the same line 
the user will enter the locations for x,y,z or for r, theta, phi in the input 
dataset, as appropriate. The locations will be entered as a nunber between 
[1, #SERIES]. 

EXAMPLE: 0,1,2, 3 i 

I 

I 
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Usage of REDO Conmand 


This command allows one to repeat a sequence of commands in batch Jobs by 
rewinding the file containing the sequence. User specifies 'REDO'. 
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Usage of the ROTATE command 


The operator rotate performs a coordinate rotation using data sets EIGVEC and 
EIGVEC1, which contain the direction cosines ( eigenvectors ) of the variance 
ellipsoid computed by operator EIG. i.e., the principal axis of the variance 
ellipsoid at a particular frequency (NFREQ). Ihus this operator assumes EIG 
hasalready been executed and the data sets EIGVEC and EIGVEC1 exist, unless 
the option to enter the direction cosines interactively is selected, in which 
case the direction cosines are entered from the terminal. 

If EIGVEC and EIGVEC 1 are used for rotation then the first three series of the 
EIGVEC dataset at NFREQ supply the direction of minimum variance. The first 
three series of the EIGVEC1 dataset at NFREQ supply the intermediate direction 
and series 4-6 of the EIGVEC 1 dataset at NFREQ supply the maximum direction of 
variation . 

If the rotation is to be performed from terminal input, the user will be 
prompted for the direction cosines. NFREQ is not required and thus the user 
will not be prompted for NFREQ. 

The name of the data set created by this operator is ROTATE.DAT. 

TO INVOKE OPERATOR; 


ROTATE datasetname 


FOR EXAMPLE: 

ROTATE IS 1201 05 THEN FOLLOW PROMPTS. 

THE PROMPTS WILL BE AS FOLLOWS: 

The first input will be to indicate where in the input dataset to find x,y,z 
(each should be a number between 1 and the number of series). The next prompt 
will indicate whether the direction cosines will be input from terminal or 
obtained from datasets EIGVEC and EIGVEC 1. If the direction cosines will be 
obtained ffom EIGVEC and EIGVEC 1 the final prompt will be for the desired 
frequency (NFREQ). If the direction cosines will be input from the terminal, 
there will be three prompts to obtain the direction cosines, 

EXAMPLE: ( DIRECTION COSINES INPUT FROM TERMINAL) 

1.2.3 1 .9890,0.2, .1439 .5567, .23, .89 . 1656, -.2367, .3456 
EXAMPLE: ( DIRECTION COSINES FROM EIGVEC AND EIGVEC1 ) 

1.2.3 0 
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Usage of SETUP Command 


SETUP receives as parameters the endpoints of the span of time that the user 
wants analyzed (IT1,1T2), and the interval size (SIZE) he wants processed* 
specified in units according to lUNIT. (Interval: basic time segment to be 
analyzed. Span: Total time composed of contiguous intervals, see Fig. 1.0-3*) 
ITl and IT2 are converted from date/time format into REAL*8 variables BEtJTIM 
and ENDTIM in units of seconds. 

At the end of the routine BEGTIM is right adjusted by the sanount of processing 
that has occurred to prepare for a subsequent run. If the left endpoint is 
greater than the right, attempt will be made to input parameters for a new 
span. 

User data are provided by calls to subroutine INPUT. Data gaps are filled with 
flag value BAD (-99*9). Analysis is based on the time sampling rate (DELTIM) 
that was found when the buffer was filled wlbh user data. 

Variables : 

IT1(6) s start date/time of span desired [4 digit year; day (Jan 1 s day 1); 
hour; min; sec; msec]. 

IT2(6) = end date/time of span desired. 

lUNIT = parameter indicating units of SIZE: (1 = seconds, 2 s minutes, 3 = 
hours, 4 = days). 

KSER = n — user specification of which series desired from INPUT block: If n is 
positive, choose the first n series in block; if n is negative, choose the 
indices represented by each digit, starting at the rightmost digit. For 
example, n = 3 chooses series 1, 2, 3; n s -351 chooses series 1, 5, 3 in that 
order. For factors influencing choice of specification order, see SPECT 
(operator which calculates spectral matrices). 

SIZE s size of interval to be processed (see lUNIT parameter). 

Those records must be placed into F0R051.DAT before processing begins, so that 
SETUP can access required parameters: 

ITl (6) format(6I5) 

KSER, SIZE, lUNIT formatd 10,F10, 3. 1 10) 

IT2(6) format(6I5). 

SETUP attempts to collect data for analysis beginning with the time specified 
as the beginning of the span. The first such point that contains all good data 
defines the start time of the interval being setup. The stop time of the' «pan 
(1T2) is an absolute maximun. The endpoint of the interval is then defined as: 

B = MIN(A4SIZE«funotion(IUNIT),IT2). 
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All points found with corresponding times less than B will be included in the 
current interval by SETUP. Ihe point corresponding to B will go into the next 
interval. Any interval may be shortened to insure that good data exists at the 
endpoint of the interval. 

For additional details on assignments and files required, see Section H and 
Appendix A and Appendix F. 

Additional span(s) subsequent to the first may be accessed by immediately 
adding the appropriate three records for each as specified above. 

Ordinarily the user need not concern himself with the size of arrays required. 
In response to the command "SETUP", the system sets up the next interval for 
processing. If the user also wishes to see a display of the intermediate data 
buffer as it is being filled, he should specify: 

SETUP 1 

If the user has his parameters stored in a file different from FOR051.DAT, he 
can access them by specifying the optional paramter "DSNsfile", for example: 

SETUP DSNsV2PUSMA.DAT 

NOTE: Do not Control Y out of IDSP if processing multiple intervals. Ihis will 
cause all of the data statements to be reinitialized and IDSP will start the 
span at the beginning. 


6.25-<2 


SETUP 4/16/84 



Usage of SHCW Coimand 


Output for every IDSP operation is placed Into a dataset identified as 
•operation.DAT'. History of any given dataset can be obtained by typxng 'SHOW 
operation'; History and data by 'SHOW operation D*. Interactive users can 
obtain hard copy of any dataset by routing output to disk via the FLOP 
command 9 and then executing 'SHOW operation D' . "D" can be any nonblank 
character . 
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Usage of SPEC! Command 


SPEC! uses the output of the Fourier Transform to compute the spectral density 
matrix (PSD). This operator uses as input the first JSER series of the 
requested dataset where JSER = min(NSER,4), and NSER is the number of series 
in the dataset. All other series are ignored. 

Let the ith point of each of the four series used as input to this operator be 
represented by A(i,J), Then PSD(i,j,k) = XNORM»Conjg( A(i J))»A(i ,k) 

gives the raw j-kth element of the ith spectral matrix. For an input dataset 
containing NPTS and for i = 1, 2, . . . , NPTS-1 XNORM r 2./( (2»NPTS-2)»»2) , 
thus folding the negative power and adding it to the positive. For i = 0, 
XNORM = 1./((2»NPTS-2)»*2). 

The user must specify how many of these matrices he wishes to have added 
together to make spectral estimates. He is required to specify an odd integer, 
thus centering the effective frequency of the estimate over the center point 
used in estimation. To retain the raw spectral estimates, the user must 
specify one point per estimate. Each estimate is actually a simple sun. 

Suppose the user specified n points per estimate where n=2k+1. The special 
case of zero frequency is accomplished by retaining the 0«th point without 
modification, and ignoring the next k points* 

If the nunber of spectral matrices was such that the final estimate is not 
"full”, the available points are used and a message is written to the history 
of the dataset, detailing how many points were used. 

By default, the system delivers results in units/HERTZ by dividing by the 
frequency width of the estimate. The 0-th estimate is skipped (left 
unchanged). If the user wishes his results in terms of power, he should 
specify the optional parameter POWER. 

The following example will form spectral matrices using dataset FFT as input 
with 7 elements forming each estimate: 

SPECT FFT 7 

The diagonal tenns 11,22,33f... are stored in SPECTD. Tiie off diagonal terms 
are stored in SPECTOFF by taking in order the rows of the upper triangular 
portion of the spectral matrix. The history record details which matrix 
element is associated with which series in dataset SPECTCFF. 


In a similar fashion, the coherence and phase arising from the off-diagonal 
terms are stored as the real and imaginary parts of complex dataset COHPH, 
respectively. WARNING: Because of the underlying mathematics, if only one 
point is used to form each estimate, the coherence will always be equal to one. 


LIMITATION— see FFT. 
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Usage of the STOP Command 


The STOP command is necessary to terminate the analysis session. A check will 
be made to insure that any subprocesses initiated for graphing have finished. 
If necessary, IDSP will wait until they finish. 


6.28-1 


STOP 11/16/84 



m 


c 



!■ >ii 


^ ' 


,''7 


Usage of SUBSER Command 

This command allows the user to form a subset of the specified dataset. The 
form of the command is: 


SUBSER dataset n1 n2 

where n1 is the index of the first series desired, and n2 is the number of 
series to be placed into the subset. 

If the dataset WINDOW has 7 series, then the following command will place 
series 2, 3* and 4 into dataset SUBSER: 

SUBSER WINDOW 2 3 

For an inverse operator, see CONCAT. 


/ 
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Usage of SUBSET Command 


The user can extract a subset of NP points from a specified dataset beginning 
at index II where II is in the interval [0;NPTS-1] and NPTS is the nunber of 
points in the parent dataset. This capability allows one to obtain contiguous 
intervals of data that has been interpolated and filtered. 

The following example extracts contiguous subsets of 500 points each from the 
dataset FILTER which has 1500 points or more. Note that the first argixnent is 
the index of the starting point, and the second argument is the nunber of 
points to be extracted. In the example we are extracting 500 points esoh time. 

SETUP 

INTERP SETUP 
FILTER INTERP 
SUBSET FILTER 0 500 

do chosen ope.ations on SUBSET(Note that 
unless COPY has been used the previous SUBSET 
has been lost) 

SL®SET FILTER 500 500 
do chosen operations on SUBSET 

SUBSET FILTER 1000 500 

If you wish to save these subsets, see COPY. 
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Usa^e of TRACE Coranand 


This command allovfs the user to form a trace of the specified dataset. The 
form of the command is: 


TRACE dataset n1 n2 

where n1 is the index of the first series desired, and n2 is the nunber of 
series to be placeo into the trace. If the dataset WINDOW has 7 series, then 
the following command will sum series 2, 3» and 4 and place the result into 
dataset TRACE: 


TRACE SPECTD 2 3 

Normally operator TRACE is used on data set SPECTD 
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Usage of WINDCW Command 


WINDOW windows every series in the dataset according to user specifications. 
(0 = Rectangular, i = 13% cosine taper, 2 = Hanning, 3 : Hamming.) 

The dataset can be padded with zeroes out to point number "n” (after 
windowing) by specifying the optional parameter "PADsn". For example, a 
Hanning window ol dataset INTERP is accomplished by: 

WINDOW INTERP 2 

A rectangular window of dataset FILDES, then padded with zeroes out to point 
2048 is accomplished by: 


WINDOW FILDES 0 PAD=2Q48 
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Usage of VTTS 

LOW PftSS FILTER WEIGHTS ESTIMATION OPERATOR 

NTS Calculates an estimate for the number of filter weights necessary to 
design a given filter using the Remez exchange algorithm. To execute, input 
the start and stop frequencies of the transition band, sampling frequency of 
the data, and the desired pass band ripple and stop band ripple. 

SAMPLE INPUTS; 

First record: F1,F2,FS 
Second record: PBR,SBR 

FI - Start frequency of the transition band 
F2 - Stop frequency of the transition band 
FS - Sampling frequency of the data 

PBR - Pass band ripple (20*log(1-«-delta. )) 

SBR - Stop band ripple (20*log(delta2)) 

Two types of ouput can be generated by WTS. First is a printout that can be 
sent to the terminal screen or to a file named [usrid.idsp3idspout.dat. To 
sent the output to a file the user must execute the command FLOP just before 
the execution of this operator. This output contains the input parameters, 
the normalized band edges, delta., delta^, the ratio between delta, and 
delta2, and the number of filter weights calculated. 

SAMPLE OUTPUT: 

TRANSITION BAND FREQUENCIES : 0.6 - 0.8 
SAMPLING FREQUENCY s 4.4 

PASS RIPPLE = 1.0 DB 
STOP RIPPLE = 50.0 DB 

BAND EDGES s 0.0000 0.1364 0.1818 0.5000 

DELTA, s 0.12202 

delta' = 0.00316 

DELTA^/DELTA2 = 38.58554 
NUMBER OF FILTER WEIGHTS = 32.459 

The second output contains the necessary input parameters for the FILDES 
command. These values will be stored in the current version of F0R058.DAT. 

SAMPLE OUTPUT: 

33,1,2,0 

0.0,0.1361,0.181 8,0.5 
1.0, 0.0 
1.0,38.586 

LIMITION: WTS will produce a FILDES input dataset with an odd number of filter 
weights (maximum of 511 weights). 


6.33-1 
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APPENDIX A Writing INPUT Routines 

The user interfaces his particular data to the IDSP system via an input 
routine which delivers a block of data to the system whenever it is needed. 

The routine must be specified as SUBROUTINE INPUT with argument list: 

( BLOCK, TIME , MAXBLK , MEMBER , NSER , DELTIM , PCNT , BAD , EOF ) . 

The routine arguments which subroutine INPUT muse supply are as follows: 

BLOCK = two-dimensional real*4 array containing up to 256 points in each of up 
to 8 series. This array must have dimension (256,8'. 

TIME = REAL*8 array containing the time associated with each point, specified 
in seconds since the beginning of I960. This array muse have dimension (256). 
To convert calandar (Jan 1 = day 1) date-time format into seconds since I960, 
a special subroutine has been supplied. Fill real*8 variable XTIME with the 
converted time by coding: 

CALL DATTIM (XTIME, YEAR, DAY, HOUR, MIN, SEC, MSEC) 

where all variables except XTIME are integer*2. YEAR is specified as a 
four-digit integer. To convert decimal days (Jan 1 r day 0) to seconds, use a 
supplied subroutine: 

call DECTIM (XTIME, year, DDAY) 

where real*8 XTIME is returned from inputs integer*2 YEAR (4 digits) and 
real*8 DDAY(Decimal Day. Fraction of Day). 

MEMBER (INTEGER*4) = number of points INPUT is delivering on this particular 
call 

NSER(INTEGER*4) = number of series being delivered to the system (max = 8) 

DELTIM (REAL*4) = time Interval in seconds between adjacent points 

PCNT(REAL*4) = percent variation user is willing to allow in DELTIM (enter It 
as .01) 

EOF s logical variable which the INPUT routine must set to .TRUE, whenever an 
end of file occurs during read, or when the user desires to inform IDSP not to 
continue calling the INPUT routine. 

The following are system arguments and must not be altered: 

MAX6LK(INTEGER*4) s maximum number of points allowed in each series s 256 

BAD(REAL*4) s IDSP flag for bad data. Any bad data detected by the input 
routine should be set to BAD (e. g. BL0CK(25,2) s BAD), Missing data should 
also be set to BAD. 

In addition, the first call to INPUT should cause character variable HISTIN to 
be initialized to a string of up to 80 characters describing the input 
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dataset. This variable must be CHARACTER*80 and be specified: 
COMMON /HIST1/ HISTIN 


For examples showing dimensions and variable declarations, see 
SYS^IDSP: INPUT. FOR, INPUTVOY.FOR, which appear in this APPENDIX. 

If the user desires to test simulated data, an easy interface is accomplished 
by writing an input routine similar to INPUTGEN.FOR, located in DRC3:[I0SP], 
replacing the statements that do the actual data generation. The user could 
make a copy of this routine and modify it to suit his needs. 

c Subroutine INPUT (stored in INPUT. FOR) 
c IDSP 

C R. M. Wenger (CSC) 3/27/81 
c 

c reads data created by GENDAT.FOR 
c Used for testing purposes. 

c 

SUBROUTINE INPUT ( BLOCK , TIME , MAXBLK , MEMBER , NSER , DELTIM , PCNT , BAD, 

1 EOF) 

COMMON /HIST1/ HISTIN 
CHARACTER»80 HISTIN 
REAL«i» BLOCK (256, 8) 

REAL»8 TIME (256), START 
LOGICAL INIT.EOF 
DATA INIT/ .FALSE./ 

IF(INIT) GO TO 75 

OPEN ( UNITS 1 1 , TYPES 'OLD' ,NAMEs ' [ZBRMW. IDSP]FOR01 1 . DAT' , 

1 FORMS 'UNFORMATTED', READONLY) 

READdDHISTIN 
READ(II) DELTIM 
INI? s .TRUE. 

EOF s .FALSE. 

PCNT s ,005 
PI s 3.1415926 
NSER s 5 
C 

75 CONTINUE 
MEMBER s 0 
DO 100 Is1,8 

READ ( 1 1 , ENDs990 ) ( BLOCK ( I , J ) , J s 1 , NSER ) , TIME ( I ) 

MEMBER s MEMBER * 1 
100 CONTINUE 
RETURN 

990 EOF s .TRUE. 

RETURN 

END 
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c — 

c Subroutine INPUTV1 (stored in INPUTVOY.FOR) 
c IDSP 
c 

c R. M. Wenger (CSC) 11/12/80. 
c 

c IMPUTV provides input from Voyager 1 summary data, 

c structure ; 

c set PCNT according to spacecraft timing variations 
c set MEMBER = 0 

c on first execution, decide which data by *'ecord of F0R052.DAT 
c loop 1 ; 

c if member + number of type .gt. MAXBLK goto Done: 
c read record 

c if EOF, skip 'extract', set EOF to true. return 

c get begin time of block in decimal day units 

c depending on type of average desired, 

call 'extract' as a function of the desired array 

increment MEMBER 
c go to loop 1 : 



c 

SUBROUTINE IN PUT ( BLOCK , TIMBLK , MAXBLK , MEMBER , NSER I , DELTIM , PCNT , BAD , 
1 EOF) 

COMMON /HIST1/ HISTIN 
CHARACTER*80 HISTIN 
REAL«4 BL0CK(256,8) 

REAL«8 TIMBLK (256) 

INTEGER*4 MEM(3) 

LOGICAL EOF.INIT 
DATA INIT /.FALSE./ 

DATA MEM/ 1, 5, 25/ 

C 

INCLUDE 'U2LJMV0Y. FOR/LIST' 

C 

IF(INIT) GO TO 100 

OPEN (UNITS 1 1 .TYPES 'OLD' .READONLY, 

1 ACCESSs 'SEQUENTIAL ' .RECORDTYPEs 'VARIABLE ' , 

2 FORMS 'UNFORMATTED' ,BL0CK3IZEs892.REC0RDSIZEs496) 

READ(52,12) HISTIN 

READ(52,11) lAV 
PCNT s .005 
EOF s .FALSE. 

INIT s .TRUE. 

100 CONTINUE 
C 

MEMBER s 0 
200 CONTINUE 

IF(MEMBER + MEM(IAV) .GT. MAXBLK) RETURN 
READ(11,ENDs600.ERRs650)(R4D(I),Is1,496) 

IFdAV .EQ. 1) CALL EXTR (H48, BLOCK, TIMBLK, IDATE.NSERI, MEM (lAV) , 

1 lAV, DELTIM, MEMBER, BAD) 

IFdAV .EQ. 2) CALL EXTR(M96.BLOCK,TIMBLK,IDATE.NSSRI,MEM(IAV) . 

1 IAV,DELTIM, MEMBER, BAD) 
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IFdAV .EQ. 3) CALL EXTR(M192,BLOCK,TIMBLK,IDATE,NSERI,MEM(IAV), 

1 lAV.DELTIM, MEMBER, BAD) 

MEMBER s MEMBER + MEM(IAV) 

GO TO 200 
600 EOF = .TRUE. 

RETURN 

650 WRITE(6,61) 
lERR = 3 
STOP 

11 FORMAT (15) 

12 FORMAT(A) 

61 F0RMAT(1X,'»**»N0TICE: ERROR READING INPUT FILES. SESSION STOPS*') 
END 
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Subroutine EXTR 


c 
c 
c 

c EXTR gets data specified in array POINTS, and places it into 

c BLOCK; it places associated times into TIMBLK. NPOINT is number 
c of points in passed array POINTS, 
c inside extract: 

c convert flag data to BAD 

c place data into BLOCK 

c calculate times (as fct of type) and place into TIMBLK 

c deliver DELTIM (in seconds) 

c — 

c 

SUBROUTINE EXTR ( POINTS , BLOCK , TIMBLK , IDATE , NSER I , NPOINT , 

1 lAV, DELTIM, MEMBER, BAD) 

REAL»4 P0INTS(275),BL0CK(256,8),DT(3) 

REAL»8 TIMBLK (256), TIME 
INTEGER»2 IDATE(6),IYR 
DATA FLAG/999./ 

DATA DT/ 48., 9.6, 1.92/ 

NSERI - 7 
DELTIM = DT(IAV) 

C 

IVR = IDATE(I) + 1900 

CALL DATTIM(TIME,IY' ,TE(2) ,IDATE(3) ,IDATE(4) ,IDATE(5) 

1 ,IDATE(o)) 

TIME = TIME - DELTIM 
DO 200 1=1, NPOINT 
TIME = TIME + DELTIM 
TIMBLKCMEMBER+I) = TIME 
NGET = I 

DC 150 Jsl.NSERI 

BLOCK(MEMBER+I,J) = POINTS(NGET) 

IF (POINTS (NGET) .KQ. FLAG) BLOCK(MEMBER+I, J) = BAD 
NGF; = NGET + NPOINT 
150 CONTINUE 
200 CONTINUE 
RETURN 
END 


V.*‘ 
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APPENDIX B-WRITING USER DEFINED OPERATORS 


The user can write operators of his own and interface them into the system. 

For how to write them, see below. Object versions of user operators can be 
interfaced into the IDSP system via the link step (see How to START Using IDSP 
Section 4.0). The sequence has been designed so that the users library is 
searched ahead of the IDSP system. 

The user can write his own operators and interface them into the IDSP system 
by adhering to the following conventions: 

1. The name of the subroutine must be "Ui" where is0,1,2,...,19. 

2. Returns must always be "RETURN 1" 

3. Use an INCLUDE statement to get current system initialization. Usage of the 
INCLUDE statement helps insure ccmipatibility with IDSP system updates. Failure 
to use INCLUDE may result in subroutine incompatibility as IDSP is updated. 

4. CMD is a character string containing user parameters specified at execution 
according to the desire of this operator. See "OPERATOR PARAMETERS VIA 
CHARACTER STRING" in APPENDIX C. 

5. HISTO is an array of character strings in which the user describes in words 
the effect of this operation on the data. 

6. HIST is the system history of operations that is retrieved with the given 
dataset when it is fetched. That history is actually a composite of all 
HISTO's that have detailed the meaning of each operation affecting the given 
dataset. Any new operation history in HISTO is appended to HIST by call to 
subroutine HISTUP. This subroutine takes NHISTO elements of array HISTO, 
appends them to the NHIST elements of array HIST and updates NHIST = NHIST + 
NHISTO. It also includes a safety check to prevent NHIST from getting too 
large. 

7. IDGET, IDSAV are variables containing the names of the datasets 
respectively that the user wishes to fetch, and then save after appropriate 
operations. 

8. Details on data structures involving subroutines FETCH and SAVE can be 
found in APPENDI’* C. 

9. For maximum utility in interactive mode, all printout except prompts should 
be written to logical unit number lOUT. This variable must not be altered by 
the subroutine. All user prompts should be written to logical unit FOR006. 

10. DEG, PI, BAD are system parameters and must not be altered and are 
respectively; 180/3.141593, 3.141593, -99.9 (BAD is flag for bad data). 

11. IDSP1 contains the space available to the user for working arrays. LIHSPC 
and MAXSER are system parameters (500,000 and 8 respectively) and must not W 
altered. 

12. The dimensions needed for arrays are calculated by subroutine LINSET and 
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are stored in array LIM. This subroutine must be called so that rrays can be 
properly dimensioned. For description of parameters, see LIMSET.FOR, APPENDIX 

C. 

The operator NORM is an example requiring one input dataset (see NORM. FOR); 
The operator EIG is an example requiring more than one input dataset (see 
EIG.FOR). 

LIST OF LOGICAL NAME ASSIGM“'=‘»"^S 
SYS$IDSP:IDSP = DRC3:[IDSP]IDSP 
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APPENDIX C-DESIGN CONSIDERATIONS FOR IDSP 


USER DATA STRUCTURE AND MANIPULATION 


1. Assume input data from spacecraft is real; usage of complex values 
occurring only after Fourit Transform. 

2. Use a single set of input/output routines to store or retrieve either real 
or complex values, with array passed as a single parameter. 

3. Each operator saves output into dataset called 'oper'.DAT. 

4. Dataset names have arbitrary length up to maximum total of thirty 
characters. 

5. User data is stored in a two-dimensional array where the first dimension 
refers to the number of points per series to facillitate operations with each 
series as a unit. 

6. Output from each operator is written unformatted to disk as data(npts,nsei ) 
regardless of whether the data are real or complex. 

7. Each operator module can set up variable domain array sizes according to 
the needs of the dataset being accessed and the working array requirements of 
the operation. The parameters NSER and NPTS are extracted from the dataset and 
used by the subroutine LIMSET to set up an array LIN defining the dimensions. 
The space available for these large arrays is found in IDSP1/SPACE. Each 
module splits this up as needed. Then these arrays and dimensions are passed 
to a routine that does the actual operation. Because of extra space 
requirements, FFT operations use a different routine LIMFFT to apportion space 
needed. For how these requirements are calculated, and limits involved, see a 
listing of the subroutine. 

8. To fapilitate implementation, each operator uses an "INCLUDE” statment: 
"INCLUDE 'OPER.FOR'" and each subprocessor uses "INCLUDE 'OPERI .FOR'". Each 
routine can add any other statments that are specific to its needs. 

9. Each operator module can specify what format it desires to have complex 

arrays in memory. If JTYPE s 1, DATA is real and must be dimensioned REAL 

DATA(LIH( 1 ) ,LIM(2) ) ; parameter JFORM is ignored. If JTYPE s 2, DATA is 

complex, and will be retrieved according to parameter JFORM. In the call to 
SUBROUTINE FETCH, if the JFORM parameter s 1, data will be returned to array 
DATA as a complex variable; DATA must be specified COMPLEX 
DATA(LIN(1) ,LIM(2>). If JFORM s 2, data will be returned to array DATA as two 
real arrays, with the real part of each series stored in DATA(npts,J) and the 
imaginary part stored in DATA(npts,nser4j) where DATA must be specified REAL 

DATA(LIN(1),LIM(2)). Similarly, data aro accessed by SUBROUTINE SAVE according 
to identification of data type and storage format given by parameters JTYPE 
and JFORM. When the calling routine is dimensioned in this fashion, the 

input/output subroutines will correctly save or fetch any of the three types 
of formats, retaining the flexibility of passing a single data array through 
an argument list. NEVER allow dimensions and NPTS and NSER to take dual roles. 
Such a structure encourages VERY SUBTLE ERRORS. 
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10. Due to special working arrays needed for operators FFT, FFTIN, and SPECT, 
a special subroutine (LIMFFT) was written to define locations and lengths of 
necessary arrays to make all these ooerators compatible. For any dataset to be 
input to any one of these three op 2 rators, it must be of size compatible with 
all three operators. 

11. All input/output of user datasets is controlled by calls to subroutines 
whose listing can be found in SYS$IDSP;I0.F0R. Further details on subroutine 
parameters are found in the comments of this li&txng. 

12. With each dataset are stored arrays .PARM and IPARM as follows: XPARM(I) = 
current data spacing XPARM(2-3) = time in seconds associated with first point 
(REAL *8). XPARM(4-5) = time in seconds associated with last point (REAL*8). 
IPARM(I) = pointer indicating current domain of dataset (1 = time (seconds); 2 
= frequency (Hertz)). IPARM(2) = version number of system (see FIXVER in 
10. FOR). 

MISCELLANEOUS DETAILS Logical File Names; 

The following FORTRAN logical unit numbers are used by IDSP and its asss/ciatfcd 
routines. 


1) User input parameters: 51-79 

LOGICAL UNIT #S 

51-79 

51 

58 

57 

2) User data: 10-19 

LOGICAL UNIT #S 
11 

13 

17 

3) System Character strings: 20-29 

LOGICAL UNIT #S 
21 
22 

23 

24 

25 


4) Error Conditions: 90-99 

LOGICAL UNIT #S 
91 


USE 

input to user INPUT routines 
input to SETUP routine 
input to FILDES (design filter) 
character strings specifying operations 


USE 

experimenter data test routine read by 

DRC3:[IDSP]INPUT.F0R 

INPUT/OUTPUT parameters 

give interactive user hardcopy 


USE 

HELP files 
command history 
submit plots 

parameters needed for separate plot 
jobs 

cluster names for plot; also used to 
associate Command Event Flag Cluster 
with a process 


USE 

error condition in Versaplot batch Job 
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Event Flags 

Event flag' 64 and 65 are reserved for testing the status of graphing jobs. 

Event flag 1 is used in GRAFCK. See operators GRAFCK, GRAPH. See programs 
GRAFGOOD, GRAFGOOF, and GRAFJB. 

IDEAS 

1. Edit wild points using IMSL routine ICSMOU 


OPERATOR PARAMETERS VIA CHARACTER STRING 


During execution, parameters for a given operation are transmitted to the 
system via a set of substrings in the character variable CMD. The command 
string for any given operator is: 

command pi p2 p3 . . . . 

where "command" is the name of the operation desired, and the "pi" are 
parameters used by tha operator. The different substrings in the variable CMD 
are always separated by blanks. Each operator is of the form "Subroutine 
opera(»,CMD)". By the time e.iecution is inside this subroutine, the "command " 
has been stripped off and CMD is left as: 

pi p2 p3 . . • 

Access to these parameters is facilitated by calling one of two subroutines: 
CALL STRING(CMD,STR1) or CALL STRNUM(CMD,NUM) . Both return the value of the 
leftmost parameter to the second argument, and retu'"n CMD with the first 
parameter and first blank stripped off, i. e., 

p<; p3 . . . 

Us6 STRING when you wish to interpret the parameter as characters; use STRNUM 
when you wish to interpret it as an integer. Usually "pi" is the name of the 
dataset desired as input for the operation. Hence, many operators CALL 
ASK(CMD,IDGET). This routine calls STRING to obtain the dataset name. If CMD 
was blank (interactive user forgot to enter dataset name before carriage 
return), the user is prompted to enter the dataset name. If the user enters an 
invalid value for an optional parameter, the system will do one of two things: 
11 issue prompt for new entry, 2) ignore the entry and use the default. 
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APPENDIX 0-Executlon of Batch Jobs 

Batch jobs require everything that interactive jobs require plus a command 
file called IDSPBAT.COM which contains the commands desired. This file must 
exist in the same directory that contains the other user parameter files. 

1) Create, for example, the following file called IDSPBAT.COM 

$ CREATE F0R057.DAT 
SETUP 

INTER P SETUP 
SHOW INTERP D 
FFT INTERP 
SHOW FFT D 
STOP 
$ EXIT 

2) On the terminal execute the command file IDSPBAT.COM 
as follows : 


eiDSPBAT.COM 

This will create the file FOR057.DAT which will be read by the batch job. 

3) Create a command file, for example BATCH.COM, containing the following: 

$ASSIGN your data FOR Oil 

any other assign statements your job may require 
$eSYS$IDSP:IDSP linknames or NOLINK 

4) Submit the coitriand file, for example BATCH.COM 


I 

i 

I 

I 

J 

) 

I 
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APPENDIX E-Conventions on DATASET NAMES 

Vfhen a dataset must be specified within a command, it can be entered in simple 
form, or can be fully qualified. Unless explicitly specified otherwise, data 
type ".DAT" is assumed. Simple form defaults to directory [usrid.IDSP]; for 
example datasetsINTERP actually specifies [usrid.IDSP3lNTERP.DAT as the 
dataset name. The user can access or copy into other directories by specifying 
something such as [usr id. ALPHA 3SETUP45 or [usrid.BETA]FFT13.SAV as the dataset 
name. The limit on the total number of characters for a dataset name is 30. 
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APPENDIX F-Exlstlng input routines 


Each application accesses data by an input routine specific to the 
experimenter data of the application. This routine is called by the system, 
and is the interface between the system and the experimenter data. Thus, user 
data is accessed by calls to subroutine INPUT which provides a block of up to 
256 points for each of up to 8 series. With each point is also delivered an 
associated time, the sane for all 8 series. Examples of existing input 
routines with required logical assignments identified by "$ASSIGN": 

NOTE; The formats of many of the tapes referred to in this Appendix can be 
found in Wish, W. H.« Commonly used digital tape, disk and card formats, 
NASA/GSFC X-694-76-242 , Re v ised February 1981. 

1) INPUTOET.FOR (reads detail data from Voyager) Prior to execution of IDSP 

the data has to be read onto disk from a Voyager Detail tape. This disk file 

is then read by INPUTDET.FOR. 

$ASSIGN experimenterdata FOR011 

Magnetic field components BX, BY, BZ are returned as series 1, 2, 3« 

Parameters must be stored in FOR052. DAT as indicated below. First record: 

string of up to 80 characters describing input file. Second record: year of 
data being accessed (4 digits) 

2) INPUT PLS. FOR (reads plasma data iVom special disk dataset prepared by 
writing only needed items— see RAYPLS.FOR) 

$ASSIGN experimenterdata FOR 011 

Eight series refer respectively to plasma data words 96-102, word 38. See p 
1.15-15 of COMMONLY USED DIGITAL TAPE, DISK AND CARD FORMATS REV FEB 1981 for 
further details. Parameters must be stored in FOR052.DAT as Indicated below. 
First record: string of up to 80 characters describing input file. 

3) INPUTVOY. FOR (reads data in Voyager Summary Tape format) 
lASSIGN experimenterdata FOR011 

Parameters must be stored in F0R052.DAT as indicated below. First record: 
string of up to 80 characters describing input file. Second record: (Specify 
in 15 type of averaged data desired: 1 s 48 sec, 2 s 9.6 sec, 3 ~ 1>92 sec) 

4) ISEEINPUT.FOR: THIS INPUT ROUTINE CAN BE USED TO READ THE ISEE-1 OR ISEE-3 
DATA POOL TAPE AND OBTAIN PARAMETERS AS REQUESTED BY THE USER FOR USE WITH THE 
IDSP PACKAGE. THE VARIABLES, FREQUENCY OF OCCURRENCE, NUMBER OF SERIES, NUMBER 
OF INTERVALS DESIRED ARE ENTERED ON A FILE F0R054.DAT. THE USER ENTERS THE 
VARIABLES DESIRED BY ENTERING THE WORD NUMBERS AS GIVEN IN "COMMONLY USED 
DIGITAL TAPE, DISK AND CARD FORMATS" PAGES 1.16-10 TO 1.16-13 FOR ISEE-3, AND 
PAGES 1.13-20 TO 1.13-27 FOR ISEE-3. FORTRAN LOGICAL UNIT 12 WILL BE USED TO 
READ THE ISEE DATA POOL INFORMATION. SEE DRA1: [U2FW0. IDSP] FOR A SAMPLE 
F0R054.DAT FILE. 
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5) RUSLIMPUT.FOR; THIS INPUT ROUTINE CAN BE USED TO READ THE ISEE-1 AND ISEE-2 

RUSSELL DATA WHICH RESIDES ON DISX AS DRC 3: [U2DHFjRUSLISEE1.DAT ISEE-1 

RUSSELL MAGNETIC FIELD DATA, AND DRC3: [U2DHF]RUSLISEE2. DAT ISEE-2 RUSSELL 

MAGNETIC FIELD DATA. THE USER ENTERS THE VARIABLES DESIRED, NUMBER CF SERIES, 
AND NUMBER OF INTERVALS DESIRED ON A FILE F0R068.DAT. THE USER ENTERS THE 
VARIABLE NUMBERS AS FOLLOWS: 09 IF WANT MAG FIELD BX COMPONENT IN MILLIGAMMAS 

(GSM), 10 IF WANT HAG FIELD BY COMPONENT IN MILLIGAMMAS (GSM), 11 IF WANT HAG 

FIELD BZ COMPONENT IN MILLIGAMMAS (GSM), 12 IF WANT NAG FIELD MAGNITUDE BT IN 
MILLIGAMMAS. FORTRAN LOGICAL UNIT 12 WILL BE USED TO READ THE RUSSELL DATA. 
SEE DRA1:CU2FW0.IDSP] FOR A SAMPLE FOR068.DAT FILE. 

6) HELIINPUT.FOR: THIS INPUT ROUTINE CAN BE USED TO READ THE HELIOS 1 AND 2 
HOURLY AVERAGE TAPE RECEIVED FROM LARRY KLEIN. THE USER ENTERS THE VARIABLE 
NUMBERS, NUMBER OF SERIES, AND THE NUMBER OF INTERVALS ON A FILE CALLED 
F0R070.DAT USING A FORMAT DESCRIPTION WRITE-UP SUPPLIED BY LARRY KLEIN. A 
SAMPLE FILE FOR070.DAT EXISTS IN DRA1: [U2FVO.I06P ]. EXAMPLE: 12 ENTERED AS 
VARIABLE NUMBER INDICATES DENSITY DESIRED. FILE FOR070.DAT WILL ALSO CONTAIN 

INFORMATION TO BE SUPPLIED ABOUT THE INPUT TAPE I.E. TAPE DRIVE, FORMAT 

TYPE, BLOCKSIZE, RECORDSIZE, FILE START, FILE STOP, AND WHETHER THE TAPE IS 

LABELLED OR NON-LABELLED. 

7) BYRNESIN.FOR: THIS INPUT ROUTINE CAN BE USED TO READ THE DATA FROM A TAPE 

SUPPLIED BY JIM BYRNES CONSISTING OF DE MAGNETIC FIELD DATA. A FILE CALLED 

F0R060.MT (SEE DRA 1: [U2FW0. IDSP] FOR AN EXAMPLE) WILL BE USED TO ENTER THE 
TAPE DRIVE, RECORD FORMAT, WHETHER LABELLED OR NOT, BLOCKSIZE, RECORDSIZE, 
START FILE, AMD STOP FILE. THIS FILE WILL ALSO CONTAIN A SECOND RECORD 

INDICATING THE NUMBER OF INTERVALS DESIRED. 

8) IMFIMPUT.FOR: THIS INPUT ROUTINE CAM BE USED TO READ THE ISEE-3 1 HOUR 

AVERAGE IMF TAPE RECEIVED FROM JOE KING. THE TAPE DRIVE USED, RECORD FORMAT, 
WHETHER SL OR NL, BLOCKSIZE, RECORDSIZE, START FILE, STOP FILE, WILL APPEAR ON 
THE FIRST RECORD OF FILE F0R061.DAT. THE SECOND RECORD WILL CONSIST OF THE 
VARIABLE NUMBERS, NUMBER CF SERIES, AND THE NUMBER OF INTERVALS. THE VARIABLE 
NUMBERS CAN BE OBTAINED FROM A FORMAT DESCRIPTION OF THE TAPE. EXAMPLE: 08 
INDICATES THAT EX IN GSE IS DESIRED. A SAMPLE F0R061.DAT FILE CAN BE SEEN IN 
DRA1: [U2FWO.I06P]. 

9) IMP8IMPUT.F0R: THIS INPUT ROUTINE CAN BE USED TO READ THE IMP H or J 15.36 

SECOND SUMMARY TAPE. VARIABLES DESIRED, NUMBER OF SERIES, NUMBER OF INTERVALS 
WILL BE ENTERED IN A FILE CALLED F0R053.DAT. THIS FILE WILL ALSO CONTAIN 

INFORMATION SUCH AS TAPE DRIVE, RECORD FORMAT, SL OR NL, BLOCKSIZE, 
RECORDSIZE, START FILE, STOP FILE. TO OBTAIN THE VARIABLE NUMBER, THE DOCUMENT 
“COMMONLY USED DIGITAL TAPE, DISK AND CARD FORMATS" PAGES 1.7-13 TO 1.7-16 
SHOULD BE CONSULTED. EXAMPLE: 10 ENTERED AS A VARIABLE NUMBER INDICATES THAT 
FIELD MAGNITUDE (FI) IS DESIRED. AN EXAMPLE OF AN F0R053.DAT FILE CAN BE SEEN 
IN DRA1: CU2FW0.IDSP]. 
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10) CMNIIMPUT.FOR: THIS INPUT ROUTINE CAN BE USED TO READ THE COMPOSITE 
OMNITAPE INFORMATION WHICH EXISTS ON DRC 3: [OMNITAPE 3 AND OBTAIN PARAMETERS AS 
REQUESTED BY THE USER FOR USE WITH THE IDSP PACKAGE. THE VARIABLES, FREQUENCY 
OF OCCURRENCE, NUMBER OF SERIES, NUMBER OF INTERVALS DESIRED ARE ENTERED ON A 
FILE F0R056.DAT. THE USER ENTERS THE VARIABLES DESIRED BY ENTERING THE WORD 
NIWBERS AS GIVEN IN "COMMONLY USED DIGITAL TAPE, DISK AND CARD FORMATS" PAGES 
1.14-1 TO 1.14-3. FORTRAN LOGICAL UNIT 12 WILL BE USED TO READ THE CMNIDATA 
RECORDS. SEE DRA1: [U2FW0.IDSP] FOR A SAMPLE F0R056.DAT FILE. 

11) VOYAINPUT.FOR; THIS INPUT ROUTINE CAN BE USED TO READ THE VOYAGER 1 AND 2 
HOURLY AVERAGE TAPE RECEIVED FROM LARRY KLEIN. THE USER ENTERS THE VARIABLE 
NUMBERS, NUMBER OF SERIES, AND THE NUMBER OF INTERVALS ON A FILE CALLED 
F0R069.DAT USING A FORMAT DESCRIPTION WRITE-UP SUPPLIED BY LARRY KLEIN. A 
SAMPLE FILE FOR069.DAT EXISTS IN DRA1: [U2FW0. IDSP]. EXAMPLE: 12 ENTERED AS 
VARIABLE NUMBER INDICATES UNSITY DESIRED. FILE F0R069.DAT WILL ALSO CONTAIN 

INFORMATION TO BE SUPPLIED ABOUT THE INPUT TAPE I.E. TAPE DRIVE, FORMAT 

TYPE, BLOCKSIZE, RECORDSIZE, FILE START, FILE STOP, AND WHETHER THE TAPE IS 
LABELLED OR NON-LABELLED. 

12) ISEEIMPIN.FOR: THIS INPUT ROUTINE CAN BE USED TO READ THE ISEE/IMP TAPE 
GENERATED BY JOE KING. THE USER ENTERS THE VARIABLES DESIRED ON A FILE 
F0R055.DAT. OlsBX OR BPAR. C2= BY OR BPER1, 3=BZ OR BPER2, 04sBM FOR ISEE-3. 
05sBX OR BPAR, 06sBY OR BPER1, OT^BZ OR BPER2, OSsEM FOR IMP. THIS FILE WILL 
ALSO CONTAIN THE NUMBER OF SERIES CHOSEN. ALSO ON A SEPARATE RECORD THIS FILE 
CONTAINS INFORMATION ABOUT THE TAPE DRIVE USED, RECORD FORMAT, SL OR NL, 
BLOCKSIZE, RECORDSIZE, START FILE, AND STOP FILE. A TAPE FORMAT DESCRIPTION 
MAY BE OBTAINED FROM JOE KING. AN EXAMPLE OF A F0R055.DAT FILE MAY BE FOUND IN 
DRA1:[U2FW0. IDSP]. 

13) ADOLVOYG.FOR: This input routine can be used to read the Voyager Conjoint 
data tape. The following parameters will be read: F2, B1, B2, B3. The user 
enters the total nimber of spans to process and whether he/she desires 1.92s, 
9.6s, or 48s data on a file called F0R067.DAT. A sample F0R)67.DAT file 
exists on DRA1: [U2FW0.IDSP]. 

14) VOYNEWINP.FOR: This input routine can be used to read the revised Voyager 
1 and 2 hourly average tape received from Larry Klein. The user enters the 
variable numbers, number of series, and the number of intervals on a file 
called FOR069.DAT using a format description write-up supplied by Larry Klein. 
A sample file F0R069.DAT exists in DRA1: [U2FW0. IDSP]. File F0R069 will also 
contain information to be supplied about the input tape, i.e., tape < ivi,, 
format type, blocksize, recordslze, file start time, file stop time, and 
whether the tape is labelled or non-labelled. 

15) YSFAHDATA.FOR: This input routine reads the file created by Fhed Herrero 
named DRA1: [YSFAH. IDSP 3NACSIN.DAT, Contact Fred Herrero for details. 

Detail on writing such a routine can be found under "Writing INPUT Routines. 
APPENDIX A". 
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APPENDIX G-Changes incorporated in Version 1 of IDSP 

The current supported version is identified by a LOGICAL assignment that is 
system generated for IDSP (see document on how to START using IDSP). The old 
version (Version 0) of IDSP will be saved but nonaccessable to the general 
user except through special arrangements with the systems management. 

New operators include EDHIST and DCL. 

Base year for internal times was changed from 1971 to I960. This time change 
was made global for the IDSP system and a special common /IDS??/ declared for 
it. The following files were affected; 

START. FOR 
SHARE. FOR 
SETUP. FOR 

IDSP was fixed to recognize datasets created under Version 0 and convert them 
internally to Version 1. All results become Version 1. This process is 
transparent to the user. A new subroutine was added to the file 10. FOR called 
FIXVER. 

The file LIMSET.FOR was modified to clarify a comment. 

The file INTERP.FOR was modified to correct a typographical error that caused 
problems for large datasets. 

The file FILTER. FOR was modified to correct a timing error that was being 
placed into the output datasets. 

The file SETUP. FOR was modified to solve several problems that had arisen. 

The operator GRAPH was modified to execute graph jobs via a subprocess rather 
than a batch job. This was done to avoid system bottleneck when the batch 
queue was. full. The STOP operator was then modified to make i.ure all graphing 
had finished before exiting IDSP. 

Appropriate changes were also made in IDSP initialization. Files affected 
were: 

GRAPH. FOR 
START. FOR 
STOP. FOR 

The file GRAFCK.FOR was modified to reflect fact that the log file of the 
graph Job is now written to the default directory. 

The file GRAFJB.FOR was modified to include the name of the dataset being 
graphed in the file GRAFJB.LOG. 

Appropriate documentation and executive modules were also modified to reflect 
these changes. 
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APPENDIX H-Existing user operators 


W. Matthaeus and M, Goldstein have written several user operators whicn may 
be of general use. Several of these are described below. As of 10/1*1/82 
these operators do not utilize the variable dimensioning aspect of IDSP. In 
most cases the data sets are dimensioned at (5002,8) although, as noted below, 
a couple of the operators use larger dimensions. The easiest way to use these 
operators is to link to DhA1:[YSWHM.IDSPU0P]USR0P.01B. 

1. BAD POINT EDITOR (operator U10) This operator allows one to scan data sets 
and remove badpoints without the necessity of using a graphics terminal. The 
program is fully interactive and the various options available are spelled out 
when one types "U10”. The program allows the user to set his own criterion 
for selecting bad points based on either the magnitude of the data or the 
standard deviation. Automatic editing based on those criteria is possible. 
In the "SCROLL" mode, four series at a time are displayed on the screen, the 
potential bad point is centered and five points before and after the selected 
point are shown. After editing the data the dataset can be rewritten on disk 
with the flagged points assigned the value "-99.9". Datasets can be edited 
rapidly with this program. 

2. SPECTRAL ANALYSIS ROUTINES (operators U11, U2, U3, U*l, U5) 

A. Spectral Analysis Using Fast Fourier Transforms (operator U11) This 
routine is similar to the SPECT operator of IDSP, but provides the 
user with more information about the statistics of the results. This 
operator is dimensioned at (20600,8); it operates on the first three 
seriec of the input dataset (usually FFT). Two output datasets are 
produced which contain the diagonal and off-diagonal elements of the 
spectral tensor as well as series which contain the magnetic 
helicity. Unlike the IDSP operator, this operator uses running 
averages to achieve whatever statistical weight the uses desires. 
Also all the information produced by the operator is written on 
dataset F0R010.DAT and can be printed after the II^P session. All 
power estimates have the dimensions of POWER in (nT) . It is assumed 
that the input data before transforming was a magnetic field time 
series in units of nT. Printed on the screen, the header of the 
output datasets, and on the F0R010.DAT dataset are the correlation 
length, the magnetic helicity length and several other length scales 
useful in MHD turbulence studies. Because this routine was written 
foi‘ analysis of solar wind data, the Taylor "frozen-in-flow" 
hypothesis is assumed and the solar wind velocity is used to convert 
from frequency to wavenumber . The output datasets contain the "raw" 
(one spectral estimate) power and magnetic helicity as well as the 
"smoothed" spectra. 

B. Blackman-Tukey "Mean Lagged Product Technique" (operators U2 U3) 
Correlation Functions (operator U2) 

The user provides a time series which, at his option, may have b.eu 
previously interpolated using the INTERP command. Fo’ lowing the 
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interactive prompts this operator produces two datasets which contain 
the diagonal and off-diagonal elements of the correlation matrix. A 
maximum lag of 10K of the total data Interval is the default 
assumption, but the user can readily request whatever he needs. 

Spectral Matrix (operator U3) The input to this operator is the 
output of operator U2. This operator WINDOWS the correlation matrix 
(Hanning is the default), then zero pads (a factor of ten is the 
default), and finally computes the fast Fourier transform. If the 
maximum lag used in operator U2 is not 10!l of the data length, then 
the amount of zero padding should be changed accordingly. One should 
know the solar wind velocity so that frequencies will be correctly 
changed to wavenumbers. The output data sets (including FOR010.DAT) 
contain the power spectrum and magnetic helicity. 

SPECTRAL ANALYSIS OF PLASMA DATA (operators U5) 

Cross helicity (operator U4) 

A quantity of great interest in MHD turbulence is the cross helicity. 
These operators permit one to determine the cross helicity using 
plasma and magnetic field data from the conjoint Voyager Summary 
tape. Operator U4 converts magnetic field measurements from units of 
nT to units of Alfven speed. The magnetic field is divided by 
SQRT(4*PI*rIio) where "rho" is the mass density. One can either use 
the instantaneous proton density on the Summary tape or the average. 
A contribution from alpha particles can be included. 

Plasma Spectra Operator (operator U5) 

U5 uses the FFT technique to calculate the power spectrum of the 
velocity, density, cross ielicity and total energy using the eight 
series that have previously been operated on by U4. Three output 
datasets are produced. 

MISC. OPERATORS 

Smoothing (operators U12 UI 3 ) Under some circumstances it is useful 
to be able to average datasets using a running average. For example, 
if one wants to perform an eigenvalue analysis one can use SPECT with 
only one spectral estimate, then smooth the spectral matrices using 
U12 (Real datasets) and U13 (Complex datasets) to produce spectra 
fully equivalent to the output of U11 but in a form that the EIG 
operator can handle. 

Operator U9 is a graph command very similar to the IDSP GRAPH but 
which allows added flexibility in choosing some of the graph 
parameters. The routine also uses smaller lettering in the headers. 
The graphs are produced interactively; no subprocess is created. 
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APPENDIX I - INSTALLATION PROCEDURE 

THE IDSP TAPE WAS CREATED USING THE DCL COPY COMMAND ON A VAX 11/780 UNDER VMS 
V2.3. THE VOLUME LABEL IS IDSPTP. THE MAGNETIC TAPE CONTAINS 124 FILES, 

DENSITYs 1600 BPI. THE TAPE CONTAINS THE FOLLOWING TYPE FILES COPIED FROM 
DISK; 

1) FORTRAN SOURCE CODE FOR IDSP ROUTINES, AND IEEE DIGITAL SIGNAL 

PROCESSING (*.FOR), 

2) COMMAND PROCEDURE FILES WHICH ARE PART OF THE IDSP PACKAGE (».COM ), 

3) DOCUMENTATION FILES WHICH ARE PART OF THE IDSP PACKAGE (».DOC), 

4) OBJECT MODULE LIBRARY OF THE IEEE DIGITAL SIGNAL PROCESSING 

(EFF , WATE, D2, GEE 2, ERROR 2, REM EZ , OUCH, 1 1MACH ) (DS P. OLB) , 

5) OBJECT MODULE LIBRARIES (F HP 2648A INTERACTIVE GRAPHICS ROUTINES 

WRITTEN BY MR. THURSTON CARLETON OF NASA/GSFC AND USED BY IDSP: 

(HPSUB.OLB AMD COMMAND. OLB) . NOTE THAT THE SOURCE CODE FOR MR. 

CARLETON 'S HP ROUTINES ARE NOT SUPPLIED, BUI' IF YOU HAVE A NEED FOR 
THEM THEY CAN BE SUPPLIED. 

6) CFJECT MODULE LIBRARY OF IDSP ROUTINES; (IDSP. OLB), 

7) OBJECT MODULE LIBRARY OF IMSL SINGLE PRECISION ROUTINES 

(EIGRS, FFTCC, UGETIO, USPKD, EHOBKS, EH0USS,EQRT2S,FFT2C, UERTST, AND 
FFTRC) aMSLIBS.OLB), 

8j two data FILES (FOR051.DAT, AND FOR060.DAT) USED BY THE SAMPLE INPUT | 

ROUTINE SUPPLIED (IDSPIN.FOR). 

9) AN OBJECT FILE OF THE IDSP. FOR ROUTINE USED BY IDSPLINK.COM AT TIME 
LINK IS DON'' THE FILE NAME IS IDSP. OBJ. 

IN SUMMARY THE TAPE HAS 6 TYPE DISK FILES COPIED ONTO IT; *.COM, ».DOC, ».FOR, 

•.OBJ, ».DAT, AND ».OLB FILES. 

I 

NOTE THAT NONE OF THE PROPRIETARY VERSAPLOT SOFTWARE HAS BEEN PROVIDED. 

THEREFORE THE IDSP OPERATORS WHICH MAKE USE OF THE VERSATEC PLOTTERS WILL NOT 
WORK. HOWEVER IF YOUR INSTALUTION HAS THE SOFTWARE AVAILABLE THEN THE 
OPERATORS SHOULD BE USEABLE. NOTE THAT COMMAND PROCEDURE GRAFJB.COM LINKS 
WHH THE VERSAPLOT SOFTWARE SO UNLESS THE SOFTWARE IS AVAILABLE THIS CWMAND 
PROCEDURE WILL NOT WORK. ALSO, COMMAND PROCEDURE GRAFCM.COM WILL NOT WORK 
UNLESS THE VERSAPLOT S(FTWARE IS AVAILABLE. IN COMMAND PROCEDURE GRAFCM.COM 
THE SYMBOL PLOT IS USED. AT OUR INSTALLATION THE SYMBOL PLOT HAS BEEN DEFINED 
AS FOLLOWS; €SYS$DISK: [SYSMCR ]PLOT WHERE THE LOGICAL NAME SYS$SYSDISK HAS BEEN 
EQUATED TO DEVICE DBAO; (PLOT EXECUTES A SYSTEM COMMAND PROCEDURE WHICH ALLOWS 
THE PLOT VECTOR FILES TO BE WRITTEN UNDER A SCRATCH UIC SO THAT USERS DO NOT 
EXCEED QUOTAS WHEN PLOTTING). 
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INSTALL IDSP THE FOLLOWING PROCEDURE COULD BE USED; 

1) CREATE A DIRECTORY (NOT SUBDIRECTORY) ON A DISK DRIVE OF YOUR CHOICE. 
FOR PURPOSE OF EXAMPLE WE ASSUME THAT A DIRECTORY [IDSPTAPE] HAS BEEN 
CREATED ON DRIVE DRA2. THIS DIRECTORY SHOULD BE ABOUT 3000 BLOCKS IN 
SIZE. 

2) EQUATE THE LOGICAL NAME SYS$IDSP TO THE DEVICE AND THE DIRECTORY 
CREATED IN ITEM #1. EXAMPLE; 

ASS DRA2; flDSPTAPE] SYS$IDSP. 

FOR THE PURPOSES OF EXAMPLE WE WILL ASSUME THAT THE LOGICAL NAME 
SYS$IDSP HAS BEEN EQUATED TO DISK DRA2 AND DIRECTORY [IDSPTAPE]. 

3) MOUNT THE TAPE ON A TAPE DRIVE; 

MOUNT MfAl; IDSPTP 

WHERE MTA1; IS THE DEVICE NAME FOR THE TAPE DRIVE AND IDSPTP IS THE 
VOLUME LABEL. A 1600 BPI COMPATIBLE TAPE DRIVE SHOULD BE USED. 

iJ) SET THE DEFAULT DIRECTORY TO SYS$IDSP. THIS ASSUMES THE DIRECTORY 
EXISTS. 

5) COPY THE TAPE ONTO SYS $I DSP AS FOLLOWS; 

COPY MTA1;».».« SYS$IDSP;«.».» 

MTA1; IS THE TAPE DRIVE. WHEN THE COPY IS COMPLETE ( SHOULD BE 12« 
FILES), 

DISMOUNT MTA1; 

(OR THE APPROPRIATE TAPE DRIVE). 

6) EDIT COMMAND PROCEDURE IDSP.COM IN SYS$IDSP AS FOLLOWS; SET UDISK 
EQUAL TO THE DISK DRIVE DEVICE NAME TO WHICH THF LOGICAL NAME 
SYS$IDSP HAS BEEN EQUATED ( LOCATED AT STATEMENT 1800 OF THE COMMAND 
PROCEDURE). AFTER MODIFYING IDSP.COM RESAVE IT. 

7) CREATE SUBDIRECTORY IDSP WITHIN THE DIRECTORY EQUATED TO BY SYS$IDSP 
AS FOLLOWS; 

CREATE/DIR DRA2; [IDSPTAPE. IDSP] 

THIS ASSUMES THAT THE DIRECTORY WAS IDSPTAPE AND THE DEVICE WAS DRA2. 

8) SET THE DEFAULT TO THE IDSP SUBDIRECTORY CREATED IN IN ITEM #7. 
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9) COPY IDSPIN.FOR, FOR051.DAT, AND F0R060.DAT INTO THE JUST CREATED 
SUBDIRECTORY ( [«.IDSP] WHERE » IS THE DIRECTORY NAME ). 

COPY SYS$IDSP; IDSPIN.FOR IDSPIN.FOR 

COPY SYS$IDGP;FOR051.DAT FOR051.DAT 

COPY SYS$IDSP:FOR060.DAT FOR06O.DAT 

FOR TEST PURPOSES, COMPILE IDSPIN.FOR TO OBTAIN IDSPIN.OBJ. IDSPIN IS 
AN "INPUT” ROUTINE THAT HAS BEEN INCLUDED TO AILOW THE JUST INSTALLED 
IDSP SYSTEM TO BE TESTED. ALSO INCLUDED ARE THE RELATED DATASETS 
F0R051.DAT AND F0R060.DAT. FOR060.DAT IS READ BY IDSPiW. F0R051.DAT 
IS READ BY THE SETUP OPERATOR ( SEE TM 83997 PAGES 6.24-142 FOR A 
DESCRIPTION OF F0R051.DAT). IDSPIN.FOR IS CURRENTLY SET UP TO 

PROVIDE PERIODIC SIMULATED DATA TO IDSP. THE CHARACTERISTICS CF THE 
DATA CAN BE CHANGED BY CHANGING FOR060.DAT AND/OR IDSPIN.FOR ITSELF. 
OF COURSE, IF YOU MODIFY IDSPIN.FOR YOU WILL HAVE TO RECOMPILE AND 
RELINK. 

10) EXECUTE THE COMMAND PROCEDURE IDSP.COM FOUND IN THE DIRECTORY EQUATED 
TO BY SYS$IDSP. SEE PACE 4.0-1 (ITEM #4) OF THE TM 83997 ON THE IDSP 
PACKAGE FOR INSTRUCTIONS ON EXECUTING THE COMMAND PROCEDURE. 
EXECUTED AS FOLLOWS; (AN EXAMPLE FOLLOWS): 

eSYS$IDSP;ID3P DRA2; [IDSPTAPE. IDSP ]IDS PIN 

11) AFTER THE IDSP SYSTEM HAS BEEN LINKED, YOU WILL SEE "ENTER COMMAND" 
ON THE SCREEN. THE FIRST COMMAND ENTERED IN THIS TEST SHOULD BE 
"SETUP" (NOTE THE IDSP COMMANDS SHOULD BE IM UPPERCASE) WHICH WILL 
READ THE TEST DATA INTO IDSP. THEN DO A "SHOW SETUP D" TO DISPLAY 
THE DATA. CONTINUE TO EXECUTE VARIOUS IDSP COMMANDS, PERHAPS THE 
COMMAND SEQUENCE SHOWN BELOW (SEE SECTION 6 OF TM 83997 FOR A 
DESCRIPTION OF THESE OPERATORS). IT IS ALSO INSTRUCTIVE TO EXECUTE 
THE IDSP "DCL" COMMAND AND THEN DO A "DIR" TO SEE THE FILES THAT HAVE 
BEEN GENERATED. 

(NOTE TO EXIT IDSP ENTER THE COMMAND STOP). 

SETUP 

SHOW SETUP D 
DSTAT SETUP 
WINDOW SETUP 0 
FFT WINDOW 
SPBCT FFT 1 

EIG SPECTD SPECTOFF DSTAT 


ETC. 


STOP 
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APPENDIX J - AN INPUT ROUTINE FOR GENERATION OF TEST DATA 


ORIGINAL PAGE iS 
OF POOR QUALITY 


SUBROUTINE INPUT(BLOCK,TIME,MAXBLK, MEMBER, NSER 
l,DELTIM,PCNT,BAD,EOF) 

COMMON/HIST 1/HISTIN 
CHARACTER»80 HISTIN 
REAL *4 BL0CK(256,8) 

REAL«8 TIME(256),XTIME,DSEED 

DATA TWOPI/6.2831853/,F1/20./,F2/18./,A1/50./,A2/1./ 
1PHI 1/0. 0/. PHI2/0. 0/.A3/1 0. / 

LOGICAL EOF, INIT/. TRUE./ 

HISTIN ='DYANMIC RANGE AND RESOLUTION TEST' 
IFdNIT.EQ. FALSE. )GOT0200 
R£AD(60, 1000)A .,A2,A3,F1,F2.F3,rai1,PHI2, PHI3 
1000 FORMAT (9F 5.0) 

MEMBER ::256 
NSER =3 

DSEEDs 123457.D0 
DELTIM=.040 
PC NT =0.001 
TsO.O 

CALL DATTIM(XTIME, 1981, 90, 0, 0, 0, 0 ) 

200 CONTINUE 

DO 100 1=1, MEMBER 
TIME(I)=T+XTIME 

BLOC K (1 , 1 )=A 1 «COS ( (TWOPI *F 1 ) »T +PH1 1 ) 

1+ A2«CaS((TWOPI»F2)«T+PHI2) 

2+ A3*COS((TWOPI*F3)«T+PHI3) 

C RAND= ( A3 *GGUBFS (DSEED ) ) 

C RAND=SIGN (RAND, BLOC K(I, D) 

C SLOCKd, ()=BLOCK(I,1)+RAND 

BLOCKd , 2)=A2«SIN ( (TWOPI *F2)«T )-A2*S’N ( (TWOPI »F3 )*T ) 
BL0CK(I,3) = 1.0 
BLOCKd, 4)=BL0CK(1, 1) 

T=T+DEI.TIM 
100 CONTINUE 

INIT=. FALSE. 

RETURN 

END 


FOR051.DAT 

1981.90.0. 0.0.0 

3 . 20 ., 1 

/J 1981,90,0,2,4,0 

F0R060.DAT 

^ 0. 10 0.20 0.40 4.70 7.80 10.9 

■I 

) 

4 


Page J-1 

APPENDIX J 3/5/84 


r 




. t ^ 





/'• / 


ABSTRACT 


The Interactive Digital Signal Processor (IDSP) is implemented on a V'X 

11/780 under VMS. It consists of a set of time series analysis "Qperatwi's" 
each of which operates on an input file to produce an output file ; the 
operators can be executed in any order that makes sense and recursively, it 
desired. The operators are the various algorithms that have been used in 
digital time series analysis work over the years. In addition, there is 

provision for user written operators to be easily Interfaced to the system. 
The system can be operated both interactively and in batch mode. 

In IDSP a file can consist of up to n (currently ns8) simultaneous time 
series. Thus storage for a file can be subdivided such that it is used, for 
example, entirely for one long single time series or for as many as n shorter 
time series, such as the components of a vector. An operator always operates 
simultaneously on all of the time series In a file. 

IDSP currently includes over thirty standard operators that range from Fourier 
transform operations ( FFT,FFTIN, WINDOW, SPECT) , design and application of 
digital filters ( FILDES,FILOPT, FILTER, WTS) , eigenvalue analysis (EIG), to 
operators that provide graphical output (GRAPH,GRAFCK) , allow batch operation 
(REDO), editing (CONCAT.EDHIST, EDIT, INTERP, SUBSET, SUBSER) and display 
information (SHOW, CNDHIS) . The complete set of standard operators is listed 
below. 

AVER, CHPHIS CONCAT, COPY, DCL, DSTAT, EDHIST, EDIT, I'in, rjLDES, FILOPT, 
FILTER, FFT, FFTIN. FLOP, GRAFCK, GRAPH, IHTERP, MF>', MNFLP, NORM, RECPOL, 
REDO, ROTATE, SETUP, SHOW, SPECT, SUUSER, SUBSET, TRACE, WINOrW, WTS, STOP. 

IDSP is being used extensively to process data sets obtained from .'cientific 
experiments onboard spacecraft such as Dynamics Explorer, ISEE, IMP and 
Voyager. In addition IDSP provides an excellent teaching tool for 

demonstrating the application of the various time series operators to 
artlfically-generated signals . 

IDSP is available from the Computer Software Management and Information Center 
(COSMIC), 112 Barrow Hall, University of Georgia, Athens, Georgia 30602, 
Program Number GSC-12862. 



